1 /*============================================================================
2 * Volume mass injection and associated source terms computation.
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 <stdio.h>
35 #include <stdlib.h>
36 #include <string.h>
37 #include <math.h>
38
39 /*----------------------------------------------------------------------------
40 * Local headers
41 *----------------------------------------------------------------------------*/
42
43 #include "bft_error.h"
44 #include "bft_mem.h"
45 #include "cs_array.h"
46 #include "cs_base.h"
47 #include "cs_equation_param.h"
48 #include "cs_field.h"
49 #include "cs_field_pointer.h"
50 #include "cs_math.h"
51 #include "cs_mesh.h"
52 #include "cs_mesh_quantities.h"
53 #include "cs_parall.h"
54 #include "cs_time_step.h"
55
56 /*----------------------------------------------------------------------------
57 * Header for the current file
58 *----------------------------------------------------------------------------*/
59
60 #include "cs_volume_mass_injection.h"
61
62 /*----------------------------------------------------------------------------
63 * Prototypes for Fortran-defined function
64 *----------------------------------------------------------------------------*/
65
66 void
67 cs_f_volume_mass_injection_get_arrays(int var_id,
68 cs_lnum_t *ncesmp,
69 cs_lnum_t **icetsm,
70 int **itpsmp,
71 cs_real_t **smcelp);
72
73 /*----------------------------------------------------------------------------*/
74
75 BEGIN_C_DECLS
76
77 /*=============================================================================
78 * Additional doxygen documentation
79 *============================================================================*/
80
81 /*!
82 \file cs_volume_mass_injection.c
83 Volume mass injection and associated source terms computation.
84
85 */
86
87 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
88
89 /*=============================================================================
90 * Local Type Definitions
91 *============================================================================*/
92
93 /*! List of cells with mass source terms */
94
95 typedef struct {
96
97 cs_lnum_t n_elts; /*!< local number of associated elements */
98 cs_lnum_t *elt_num; /*!< local cell ids, (1-based for Fortran
99 compatibilty, as it may appear in Fortran
100 use code, but not C user code) */
101
102 } cs_volume_mass_injection_t;
103
104 /*============================================================================
105 * Prototypes for functions intended for use only by Fortran wrappers.
106 * (descriptions follow, with function bodies).
107 *============================================================================*/
108
109 void
110 cs_f_mass_source_terms_get_pointers(cs_lnum_t *ncesmp,
111 cs_lnum_t **icetsm);
112
113 /*============================================================================
114 * Global variables
115 *============================================================================*/
116
117 cs_volume_mass_injection_t *_mass_injection = NULL;
118
119 /*============================================================================
120 * Private function definitions
121 *============================================================================*/
122
123 /*----------------------------------------------------------------------------*/
124 /*!
125 * \brief Evaluate a symmeytric-tensor-valued quantity for a list of elements
126 * This function complies with the generic function type defined as
127 * cs_xdef_eval_t
128 *
129 * \param[in] n_elts number of elements to consider
130 * \param[in] elt_ids list of element ids
131 * \param[in] dense_output perform an indirection for output (true/false)
132 * \param[in] mesh pointer to a cs_mesh_t structure
133 * \param[in] connect pointer to a cs_cdo_connect_t structure
134 * \param[in] quant pointer to a cs_cdo_quantities_t structure
135 * \param[in] time_eval physical time at which one evaluates the term
136 * \param[in] context NULL or pointer to a context structure
137 * \param[in, out] eval array storing the result (must be allocated)
138 */
139 /*----------------------------------------------------------------------------*/
140
141 static void
_eval_sym_tensor_by_val(cs_lnum_t n_elts,const cs_lnum_t * elt_ids,bool dense_output,void * context,cs_real_t * eval)142 _eval_sym_tensor_by_val(cs_lnum_t n_elts,
143 const cs_lnum_t *elt_ids,
144 bool dense_output,
145 void *context,
146 cs_real_t *eval)
147 {
148 if (n_elts == 0)
149 return;
150
151 const cs_real_t *constant_val = (cs_real_t *)context;
152
153 /* Sanity checks */
154 assert(eval != NULL && constant_val != NULL);
155
156 if (elt_ids != NULL && !dense_output) {
157
158 # pragma omp parallel for if (n_elts > CS_THR_MIN)
159 for (cs_lnum_t i = 0; i < n_elts; i++) {
160
161 cs_real_t *_res = eval + 6*elt_ids[i];
162
163 _res[0] = constant_val[0];
164 _res[1] = constant_val[1];
165 _res[2] = constant_val[2];
166 _res[3] = constant_val[3];
167 _res[4] = constant_val[4];
168 _res[5] = constant_val[5];
169
170 }
171
172 }
173 else {
174
175 # pragma omp parallel for if (n_elts > CS_THR_MIN)
176 for (cs_lnum_t i = 0; i < n_elts; i++) {
177
178 cs_real_t *_res = eval + 3*i;
179
180 _res[0] = constant_val[0];
181 _res[1] = constant_val[1];
182 _res[2] = constant_val[2];
183 _res[3] = constant_val[3];
184 _res[4] = constant_val[4];
185 _res[5] = constant_val[5];
186
187 }
188
189 }
190 }
191
192 /*----------------------------------------------------------------------------*/
193 /*!
194 * \brief Evaluate a symmetric-tensor-valued quantity with only a
195 * time-dependent variation for a list of elements
196 * This function complies with the generic function type defined as
197 * cs_xdef_eval_t
198 *
199 * \param[in] n_elts number of elements to consider
200 * \param[in] elt_ids list of element ids
201 * \param[in] dense_output perform an indirection for output (true/false)
202 * \param[in] mesh pointer to a cs_mesh_t structure
203 * \param[in] connect pointer to a cs_cdo_connect_t structure
204 * \param[in] quant pointer to a cs_cdo_quantities_t structure
205 * \param[in] time_eval physical time at which one evaluates the term
206 * \param[in] context NULL or pointer to a context structure
207 * \param[in, out] eval array storing the result (must be allocated)
208 */
209 /*----------------------------------------------------------------------------*/
210
211 static void
_eval_sym_tensor_at_cells_by_time_func(cs_lnum_t n_elts,const cs_lnum_t * elt_ids,bool dense_output,cs_real_t time_eval,void * context,cs_real_t * eval)212 _eval_sym_tensor_at_cells_by_time_func(cs_lnum_t n_elts,
213 const cs_lnum_t *elt_ids,
214 bool dense_output,
215 cs_real_t time_eval,
216 void *context,
217 cs_real_t *eval)
218 {
219 cs_xdef_time_func_context_t *tfc = (cs_xdef_time_func_context_t *)context;
220
221 /* Sanity checks */
222 assert(tfc != NULL);
223
224 /* Evaluate the quantity */
225 cs_real_t _eval[6];
226 tfc->func(time_eval, tfc->input, _eval);
227
228 if (elt_ids != NULL && !dense_output) {
229
230 # pragma omp parallel for if (n_elts > CS_THR_MIN)
231 for (cs_lnum_t i = 0; i < n_elts; i++)
232 for (int k = 0; k < 6; k++)
233 eval[6*elt_ids[i] + k] = _eval[k];
234
235 }
236 else {
237
238 # pragma omp parallel for if (n_elts > CS_THR_MIN)
239 for (cs_lnum_t i = 0; i < n_elts; i++)
240 for (int k = 0; k < 6; k++)
241 eval[6*i+k] = _eval[k];
242
243 }
244 }
245
246 /*----------------------------------------------------------------------------*/
247 /*!
248 * \brief Implicit and explicit mass source terms computation.
249 *
250 * \param[in,out] xdef volume injection definition
251 * \param[out] st_loc explicit source term part independant
252 * of the variable
253 */
254 /*----------------------------------------------------------------------------*/
255
256 static void
_volume_mass_injection_eval(cs_xdef_t * def,cs_real_t st_loc[])257 _volume_mass_injection_eval(cs_xdef_t *def,
258 cs_real_t st_loc[])
259 {
260 const cs_mesh_t *m = cs_glob_mesh;
261 const cs_real_t t_eval = cs_glob_time_step->t_cur;
262
263 const cs_zone_t *z = cs_volume_zone_by_id(def->z_id);
264
265 const cs_cdo_connect_t *connect = NULL;
266 const cs_cdo_quantities_t *quant = NULL;
267
268 bool dense = true; /* dense (zone-based) here */
269
270 cs_lnum_t dim = def->dim;
271
272 switch(def->type) {
273 case CS_XDEF_BY_ANALYTIC_FUNCTION:
274 cs_xdef_eval_at_cells_by_analytic(z->n_elts,
275 z->elt_ids,
276 dense,
277 m,
278 connect,
279 quant,
280 t_eval,
281 def->context,
282 st_loc);
283 break;
284
285 case CS_XDEF_BY_ARRAY:
286 cs_xdef_eval_scalar_at_cells_by_array(z->n_elts,
287 z->elt_ids,
288 dense,
289 m,
290 connect,
291 quant,
292 t_eval,
293 def->context,
294 st_loc);
295 break;
296
297 case CS_XDEF_BY_DOF_FUNCTION:
298 bft_error(__FILE__, __LINE__, 0,
299 _("In %s, evaluation by DOF function not supported"),
300 __func__);
301 break;
302
303 case CS_XDEF_BY_FIELD:
304 cs_xdef_eval_cell_by_field(z->n_elts,
305 z->elt_ids,
306 dense,
307 m,
308 connect,
309 quant,
310 t_eval,
311 def->context,
312 st_loc);
313 break;
314
315 case CS_XDEF_BY_QOV:
316 {
317 if (dim == 1) {
318 cs_xdef_eval_scalar_by_val(z->n_elts,
319 z->elt_ids,
320 dense,
321 m,
322 connect,
323 quant,
324 t_eval,
325 def->context,
326 st_loc);
327 for (cs_lnum_t i = 0; i < z->n_elts; i++)
328 st_loc[i] /= z->f_measure;
329 }
330 else if (dim == 3) {
331 cs_xdef_eval_vector_by_val(z->n_elts,
332 z->elt_ids,
333 dense,
334 m,
335 connect,
336 quant,
337 t_eval,
338 def->context,
339 st_loc);
340 for (cs_lnum_t i = 0; i < z->n_elts; i++) {
341 for (cs_lnum_t j = 0; j < 3; j++)
342 st_loc[i*3 + j] /= z->f_measure;
343 }
344 }
345 else if (dim == 6) {
346 _eval_sym_tensor_by_val(z->n_elts,
347 z->elt_ids,
348 dense,
349 def->context,
350 st_loc);
351 for (cs_lnum_t i = 0; i < z->n_elts; i++) {
352 for (cs_lnum_t j = 0; j < 6; j++)
353 st_loc[i*6 + j] /= z->f_measure;
354 }
355 }
356 }
357 break;
358
359 case CS_XDEF_BY_SUB_DEFINITIONS:
360 bft_error(__FILE__, __LINE__, 0,
361 _("In %s, evaluation by sub-definitions not supported"),
362 __func__);
363 break;
364
365 case CS_XDEF_BY_TIME_FUNCTION:
366 if (dim == 1)
367 cs_xdef_eval_scalar_at_cells_by_time_func(z->n_elts,
368 z->elt_ids,
369 dense,
370 m,
371 connect,
372 quant,
373 t_eval,
374 def->context,
375 st_loc);
376 else if (dim == 3)
377 cs_xdef_eval_vector_at_cells_by_time_func(z->n_elts,
378 z->elt_ids,
379 dense,
380 m,
381 connect,
382 quant,
383 t_eval,
384 def->context,
385 st_loc);
386 else if (dim == 6)
387 _eval_sym_tensor_at_cells_by_time_func(z->n_elts,
388 z->elt_ids,
389 dense,
390 t_eval,
391 def->context,
392 st_loc);
393 break;
394
395 case CS_XDEF_BY_VALUE:
396 if (dim == 1)
397 cs_xdef_eval_scalar_by_val(z->n_elts,
398 z->elt_ids,
399 dense,
400 m,
401 connect,
402 quant,
403 t_eval,
404 def->context,
405 st_loc);
406 else if (dim == 3)
407 cs_xdef_eval_vector_by_val(z->n_elts,
408 z->elt_ids,
409 dense,
410 m,
411 connect,
412 quant,
413 t_eval,
414 def->context,
415 st_loc);
416 else if (dim == 6)
417 _eval_sym_tensor_by_val(z->n_elts,
418 z->elt_ids,
419 dense,
420 def->context,
421 st_loc);
422 break;
423
424 default:
425 bft_error(__FILE__, __LINE__, 0,
426 _("In %s, cs_xdef_t type %d not supported"),
427 __func__, def->type);
428 }
429 }
430
431 /*============================================================================
432 * Fortran wrapper function definitions
433 *============================================================================*/
434
435 /*----------------------------------------------------------------------------*/
436 /*!
437 * \brief Get pointer to mass source terms cell count and list.
438 *
439 * \param[out] ncesmp number of cells with mass source terms
440 * \param[out] icetsm list of cells with mass source terms (1-based numbers)
441 */
442 /*----------------------------------------------------------------------------*/
443
444 void
cs_f_mass_source_terms_get_pointers(cs_lnum_t * ncesmp,cs_lnum_t ** icetsm)445 cs_f_mass_source_terms_get_pointers(cs_lnum_t *ncesmp,
446 cs_lnum_t **icetsm)
447 {
448 *ncesmp = 0;
449 *icetsm = NULL;
450
451 if (_mass_injection != NULL) {
452 *ncesmp = _mass_injection->n_elts;
453 *icetsm = _mass_injection->elt_num;
454 }
455 }
456
457 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
458
459 /*============================================================================
460 * Public function definitions
461 *============================================================================*/
462
463 /*----------------------------------------------------------------------------*/
464 /*!
465 * \brief Flag volume zones with the appropriate
466 * CS_VOLUME_ZONE_MASS_SOURCE_TERM flag when at least one volume
467 * mass injection on that zone is present.
468 *
469 * This is necessary for the reverse zone indexing required by the legacy code
470 * to function with defintions that are partially unrolled an not purely
471 * zone-based.
472 */
473 /*----------------------------------------------------------------------------*/
474
475 void
cs_volume_mass_injection_flag_zones(void)476 cs_volume_mass_injection_flag_zones(void)
477 {
478 /* First, flag all volume zones with injection definitions */
479
480 /* Now add cells in injection zones based on pressure
481 (injection values of other variables where no definition is
482 associated to the pressure may be ignored as they would be
483 multiplied by zero in any case) */
484
485 cs_field_t *f = cs_field_by_name_try("pressure");
486
487 if (f != NULL) {
488 if (! (f->type & CS_FIELD_VARIABLE))
489 f = NULL;
490 }
491
492 if (f != NULL) {
493
494 /* Retrieve the equation param to set */
495
496 cs_equation_param_t *eqp
497 = cs_field_get_key_struct_ptr(f, cs_field_key_id("var_cal_opt"));
498
499 for (int i = 0; i < eqp->n_volume_mass_injections; i++) {
500 cs_xdef_t *v_inj = eqp->volume_mass_injections[i];
501 const cs_zone_t *z = cs_volume_zone_by_id(v_inj->z_id);
502 cs_volume_zone_set_type(z->id, CS_VOLUME_ZONE_MASS_SOURCE_TERM);
503 }
504
505 }
506 }
507
508 /*----------------------------------------------------------------------------*/
509 /*!
510 * \brief Build the list and zone ids of cells with volume mass injection.
511 *
512 * \param[in] n_cells number of cells in mass source term zones
513 * \param[out] cell_num numbers (1-based) cells in mass source term zones
514 * \param[out] cell_zone_id associated zone ids
515 */
516 /*----------------------------------------------------------------------------*/
517
518 void
cs_volume_mass_injection_build_lists(cs_lnum_t n_cells,cs_lnum_t cell_num[],int cell_zone_id[])519 cs_volume_mass_injection_build_lists(cs_lnum_t n_cells,
520 cs_lnum_t cell_num[],
521 int cell_zone_id[])
522 {
523 CS_UNUSED(n_cells); /* Avoid a warning when compiling with optimization */
524
525 /* First, flag all volume zones with injection definitions */
526
527 cs_lnum_t l = 0;
528
529 for (int z_id = 0; z_id < cs_volume_zone_n_zones(); z_id++) {
530
531 const cs_zone_t *z = cs_volume_zone_by_id(z_id);
532
533 if (! (z->type & CS_VOLUME_ZONE_MASS_SOURCE_TERM))
534 continue;
535
536 for (cs_lnum_t j = 0; j < z->n_elts; j++) {
537 cell_num[l] = z->elt_ids[j] + 1;
538 cell_zone_id[l] = z_id;
539 l++;
540 }
541
542 }
543
544 assert(l == n_cells);
545 }
546
547 /*----------------------------------------------------------------------------*/
548 /*!
549 * \brief Evaluate contributions to volume mass injection.
550 *
551 * \param[in] nvar total number of variables
552 * \param[in] ncesmp number of cells with mass source term
553 * \param[in] itypsm mass source type for the working variable
554 * size: [nvar][ncesmp]
555 * \param[in] smacel values of the variables associated to the
556 * mass source (for the pressure variable,
557 * smacel is the mass flux)
558 * size: [nvar][ncesmp]
559 */
560 /*----------------------------------------------------------------------------*/
561
562 void
cs_volume_mass_injection_eval(int nvar,cs_lnum_t ncesmp,int itypsm[],cs_real_t smacel[])563 cs_volume_mass_injection_eval(int nvar,
564 cs_lnum_t ncesmp,
565 int itypsm[],
566 cs_real_t smacel[])
567 {
568 const int key_eqp_id = cs_field_key_id("var_cal_opt");
569 const int var_key_id = cs_field_key_id("variable_id");
570
571 /* Initialize arrays */
572
573 for (cs_lnum_t ivar = 0; ivar < nvar; ivar++) {
574 for (cs_lnum_t i = 0; i < ncesmp; i++) {
575 itypsm[ncesmp*ivar + i] = 0;
576 smacel[ncesmp*ivar + i] = 0.;
577 }
578 }
579
580 /* Compute shift for zones in case they do not appear in order */
581
582 int n_zones = cs_volume_zone_n_zones();
583 cs_lnum_t *z_shift;
584 BFT_MALLOC(z_shift, n_zones, cs_lnum_t);
585
586 cs_lnum_t c_shift = 0;
587
588 for (int z_id = 0; z_id < n_zones; z_id++) {
589 const cs_zone_t *z = cs_volume_zone_by_id(z_id);
590 if (z->type & CS_VOLUME_ZONE_MASS_SOURCE_TERM) {
591 z_shift[z_id] = c_shift;
592 c_shift += z->n_elts;
593 }
594 else
595 z_shift[z_id] = -1;
596 }
597
598 int n_fields = cs_field_n_fields();
599
600 for (int f_id = 0; f_id < n_fields; f_id++) {
601
602 cs_field_t *f = cs_field_by_id(f_id);
603
604 if (! (f->type & CS_FIELD_VARIABLE))
605 continue;
606
607 /* Retrieve the equation param to set */
608
609 cs_equation_param_t *eqp = cs_field_get_key_struct_ptr(f, key_eqp_id);
610 int ivar = cs_field_get_key_int(f, var_key_id) - 1;
611
612 /* xdef-based method */
613
614 for (int inj_idx = 0; inj_idx < eqp->n_volume_mass_injections; inj_idx++) {
615
616 cs_xdef_t *v_inj = eqp->volume_mass_injections[inj_idx];
617 const cs_zone_t *z = cs_volume_zone_by_id(v_inj->z_id);
618
619 c_shift = z_shift[z->id];
620
621 if (c_shift < 0)
622 continue; /* Ignore injection if it does not appear for pressure */
623
624 const cs_lnum_t n_z_vals = z->n_elts*(cs_lnum_t)(f->dim);
625
626 cs_real_t *st_loc;
627 BFT_MALLOC(st_loc, n_z_vals, cs_real_t);
628 for (cs_lnum_t j = 0; j < n_z_vals; j++)
629 st_loc[j] = 0;
630
631 _volume_mass_injection_eval(v_inj, st_loc);
632
633 if (f->dim == 1) {
634 for (cs_lnum_t i = 0; i < z->n_elts; i++) {
635 cs_lnum_t j = c_shift + i;
636 itypsm[ivar*ncesmp + j] = 1;
637 smacel[ivar*ncesmp + j] += st_loc[i];
638 }
639 }
640 else {
641 const cs_lnum_t dim = f->dim;
642 for (cs_lnum_t i = 0; i < z->n_elts; i++) {
643 cs_lnum_t j = c_shift + i;
644 itypsm[ivar*ncesmp + j] = 1;
645 for (cs_lnum_t k = 0; k < dim; k++)
646 smacel[(ivar+k)*ncesmp + j] += st_loc[i*dim + k];
647 }
648 }
649
650 BFT_FREE(st_loc);
651
652 }
653
654 }
655
656 BFT_FREE(z_shift);
657 }
658
659 /*----------------------------------------------------------------------------*/
660 /*!
661 * \brief Return pointers to the mass source term arrays.
662 *
663 * \param[in] f pointer to associated field
664 * \param[out] ncesmp number of cells with mass source terms
665 * \param[out] icetsm pointet to source mass cells list (1-based numbering)
666 * \param[out] itpsmp mass source type for the working variable
667 * (see \ref cs_user_mass_source_terms)
668 * \param[out] s_type mass source types (0: ambient value, 1: s_val value)
669 * \param[out] smcelp pointer to mass source values
670 * \param[out] gamma pointer to flow mass value
671 */
672 /*----------------------------------------------------------------------------*/
673
674 void
cs_volume_mass_injection_get_arrays(const cs_field_t * f,cs_lnum_t * ncesmp,cs_lnum_t ** icetsm,int ** itpsmp,cs_real_t ** smcelp,cs_real_t ** gamma)675 cs_volume_mass_injection_get_arrays(const cs_field_t *f,
676 cs_lnum_t *ncesmp,
677 cs_lnum_t **icetsm,
678 int **itpsmp,
679 cs_real_t **smcelp,
680 cs_real_t **gamma)
681 {
682 *ncesmp = 0;
683 *icetsm = NULL;
684 *itpsmp = NULL;
685 *smcelp = NULL;
686 *gamma = NULL;
687
688 const int k_variable_id = cs_field_key_id("variable_id");
689
690 const int var_id = cs_field_get_key_int(f, k_variable_id);
691
692 cs_f_volume_mass_injection_get_arrays(var_id,
693 ncesmp, icetsm, itpsmp, smcelp);
694
695 if (*ncesmp > 0) {
696 cs_lnum_t _ncesmp;
697 cs_lnum_t *_icetsm;
698 int *_itpsmp;
699 const int p_var_id = cs_field_get_key_int(CS_F_(p), k_variable_id);
700 cs_f_volume_mass_injection_get_arrays(p_var_id,
701 &_ncesmp, &_icetsm, &_itpsmp,
702 gamma);
703 }
704 }
705
706 /*----------------------------------------------------------------------------*/
707
708 END_C_DECLS
709