1 #ifndef __CS_BOUNDARY_CONDITIONS_H__
2 #define __CS_BOUNDARY_CONDITIONS_H__
3
4 /*============================================================================
5 * Boundary condition handling.
6 *============================================================================*/
7
8 /*
9 This file is part of Code_Saturne, a general-purpose CFD tool.
10
11 Copyright (C) 1998-2021 EDF S.A.
12
13 This program is free software; you can redistribute it and/or modify it under
14 the terms of the GNU General Public License as published by the Free Software
15 Foundation; either version 2 of the License, or (at your option) any later
16 version.
17
18 This program is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
20 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
21 details.
22
23 You should have received a copy of the GNU General Public License along with
24 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
25 Street, Fifth Floor, Boston, MA 02110-1301, USA.
26 */
27
28 /*----------------------------------------------------------------------------*/
29
30 /*----------------------------------------------------------------------------
31 * Standard C library headers
32 *----------------------------------------------------------------------------*/
33
34 /*----------------------------------------------------------------------------
35 * Local headers
36 *----------------------------------------------------------------------------*/
37
38 #include <ple_locator.h>
39
40 #include "fvm_nodal.h"
41 #include "fvm_writer.h"
42
43 #include "cs_base.h"
44 #include "cs_field.h"
45 #include "cs_math.h"
46 #include "cs_mesh_location.h"
47
48 /*----------------------------------------------------------------------------*/
49
50 BEGIN_C_DECLS
51
52 /*============================================================================
53 * Macro definitions
54 *============================================================================*/
55
56 /*============================================================================
57 * Local type definitions
58 *============================================================================*/
59
60 /*=============================================================================
61 * Global variables
62 *============================================================================*/
63
64 /*! Boundary condition type (code) associated with each boundary face */
65
66 extern const int *cs_glob_bc_type;
67
68 /*! boundary zone number associated with each boundary face
69 (specific physical models)*/
70
71 extern const int *cs_glob_bc_face_zone;
72
73 /*============================================================================
74 * Public function prototypes
75 *============================================================================*/
76
77 /*----------------------------------------------------------------------------
78 * Handling of boundary condition definition errors and associated output.
79 *
80 * This function checks for errors, and simply returns if no error is
81 * encountered. In case of error, it outputs helpful information so as to
82 * make it easier to locate the matching faces.
83 *
84 * For each boundary face, bc_type defines the boundary condition type.
85 * As a convention here, zero values correspond to undefined types,
86 * positive values to defined types (with no error), and negative values
87 * to defined types with inconsistent or incompatible values, the
88 * absolute value indicating the original boundary condition type.
89 *
90 * An optional label may be used if the error is related to another
91 * attribute than the boundary type, for appropriate error reporting.
92 *
93 * parameters:
94 * bc_flag <-- array of BC type ids
95 * type_name <-- name of attribute in error, or NULL
96 *----------------------------------------------------------------------------*/
97
98 void
99 cs_boundary_conditions_error(const int *bc_type,
100 const char *type_name);
101
102 /*----------------------------------------------------------------------------
103 * Locate shifted boundary face coordinates on possibly filtered
104 * cells or boundary faces for later interpolation.
105 *
106 * parameters:
107 * location_type <-- matching values location (CS_MESH_LOCATION_CELLS
108 * or CS_MESH_LOCATION_BOUNDARY_FACES)
109 * n_location_elts <-- number of selected location elements
110 * n_faces <-- number of selected boundary faces
111 * location_elts <-- list of selected location elements (0 to n-1),
112 * or NULL if no indirection is needed
113 * faces <-- list of selected boundary faces (0 to n-1),
114 * or NULL if no indirection is needed
115 * coord_shift <-- array of coordinates shift relative to selected
116 * boundary faces
117 * coord_stride <-- access stride in coord_shift: 0 for uniform
118 * shift, 1 for "per face" shift.
119 * tolerance <-- relative tolerance for point location.
120 *
121 * returns:
122 * associated locator structure
123 *----------------------------------------------------------------------------*/
124
125 ple_locator_t *
126 cs_boundary_conditions_map(cs_mesh_location_type_t location_type,
127 cs_lnum_t n_location_elts,
128 cs_lnum_t n_faces,
129 const cs_lnum_t *location_elts,
130 const cs_lnum_t *faces,
131 cs_real_3_t *coord_shift,
132 int coord_stride,
133 double tolerance);
134
135 /*----------------------------------------------------------------------------
136 * Set mapped boundary conditions for a given field and mapping locator.
137 *
138 * parameters:
139 * field <-- field whose boundary conditions are set
140 * locator <-- associated mapping locator, as returned
141 * by cs_boundary_conditions_map().
142 * location_type <-- matching values location (CS_MESH_LOCATION_CELLS or
143 * CS_MESH_LOCATION_BOUNDARY_FACES)
144 * normalize <-- normalization option:
145 * 0: values are simply mapped
146 * 1: values are mapped, then multiplied
147 * by a constant factor so that their
148 * surface integral on selected faces
149 * is preserved (relative to the
150 * input values)
151 * 2: as 1, but with a boundary-defined
152 * weight, defined by balance_w
153 * 3: as 1, but with a cell-defined
154 * weight, defined by balance_w
155 * interpolate <-- interpolation option:
156 * 0: values are simply based on matching
157 * cell or face center values
158 * 1: values are based on matching cell
159 * or face center values, corrected
160 * by gradient interpolation
161 * n_faces <-- number of selected boundary faces
162 * faces <-- list of selected boundary faces (0 to n-1),
163 * or NULL if no indirection is needed
164 * balance_w <-- optional balance weight, or NULL
165 * nvar <-- number of variables requiring BC's
166 * rcodcl <-> boundary condition values
167 *----------------------------------------------------------------------------*/
168
169 void
170 cs_boundary_conditions_mapped_set(const cs_field_t *f,
171 ple_locator_t *locator,
172 cs_mesh_location_type_t location_type,
173 int normalize,
174 int interpolate,
175 cs_lnum_t n_faces,
176 const cs_lnum_t *faces,
177 cs_real_t *balance_w,
178 int nvar,
179 cs_real_t rcodcl[]);
180
181 /*----------------------------------------------------------------------------*/
182 /*!
183 * \brief Add location of locate shifted boundary face coordinates on
184 * cells or boundary faces for automatic interpolation.
185 *
186 * \note
187 * This function is currently restricted to mapping of boundary face
188 * locations (usually from boundary zones) to cell of boundary face
189 * locations, but could be extended to other location types in the future.
190 *
191 * \param[in] bc_location_id id of selected boundary mesh location;
192 * currently restricted to subsets of
193 * boundary faces (i.e. boundary zone
194 * location ids).
195 * \param[in] source_location_id id of selected location mesh location
196 * (usually CS_MESH_LOCATION_CELLS but can be
197 * a more restricted cell or boundary face zone
198 * location location id).
199 * \param[in] coord_shift coordinates shift relative to selected
200 * boundary faces
201 * \param[in] tolerance relative tolerance for point location.
202 *
203 * \return id of added map
204 */
205 /*----------------------------------------------------------------------------*/
206
207 int
208 cs_boundary_conditions_add_map(int bc_location_id,
209 int source_location_id,
210 cs_real_t coord_shift[3],
211 double tolerance);
212
213 /*----------------------------------------------------------------------------
214 * Create the boundary conditions face type and face zone arrays
215 *----------------------------------------------------------------------------*/
216
217 void
218 cs_boundary_conditions_create(void);
219
220 /*----------------------------------------------------------------------------
221 * Free the boundary conditions face type and face zone arrays.
222 *
223 * This also frees boundary condition mappings which may have been defined.
224 *----------------------------------------------------------------------------*/
225
226 void
227 cs_boundary_conditions_free(void);
228
229 /*----------------------------------------------------------------------------*/
230 /*!
231 * \brief Set Neumann BC for a scalar for a given face.
232 *
233 * \param[out] a explicit BC coefficient for gradients
234 * \param[out] af explicit BC coefficient for diffusive flux
235 * \param[out] b implicit BC coefficient for gradients
236 * \param[out] bf implicit BC coefficient for diffusive flux
237 * \param[in] qimp flux value to impose
238 * \param[in] hint internal exchange coefficient
239 */
240 /*----------------------------------------------------------------------------*/
241
242 inline static void
cs_boundary_conditions_set_neumann_scalar(cs_real_t * a,cs_real_t * af,cs_real_t * b,cs_real_t * bf,cs_real_t qimp,cs_real_t hint)243 cs_boundary_conditions_set_neumann_scalar(cs_real_t *a,
244 cs_real_t *af,
245 cs_real_t *b,
246 cs_real_t *bf,
247 cs_real_t qimp,
248 cs_real_t hint)
249 {
250 /* Gradient BCs */
251 *a = -qimp/hint;
252 *b = 1.;
253
254 /* Flux BCs */
255 *af = qimp;
256 *bf = 0.;
257 }
258
259 /*----------------------------------------------------------------------------*/
260 /*!
261 * \brief Set Neumann BC for a scalar for a given face.
262 *
263 * \param[out] a explicit BC coefficient for gradients
264 * \param[out] af explicit BC coefficient for diffusive flux
265 * \param[out] b implicit BC coefficient for gradients
266 * \param[out] bf implicit BC coefficient for diffusive flux
267 * \param[in] qimpv flux value to impose
268 * \param[in] hint internal exchange coefficient
269 */
270 /*----------------------------------------------------------------------------*/
271
272 inline static void
cs_boundary_conditions_set_neumann_vector(cs_real_t a[3],cs_real_t af[3],cs_real_t b[3][3],cs_real_t bf[3][3],const cs_real_t qimpv[3],cs_real_t hint)273 cs_boundary_conditions_set_neumann_vector(cs_real_t a[3],
274 cs_real_t af[3],
275 cs_real_t b[3][3],
276 cs_real_t bf[3][3],
277 const cs_real_t qimpv[3],
278 cs_real_t hint)
279 {
280 /* Gradient BCs */
281
282 for (size_t i = 0; i < 3; i++) {
283 a[i] = -qimpv[i] / fmax(hint, 1.e-300);
284 }
285
286 b[0][0] = 1., b[0][1] = 0., b[0][2] = 0.;
287 b[1][0] = 0., b[1][1] = 1., b[1][2] = 0.;
288 b[2][0] = 0., b[2][1] = 0., b[2][2] = 1.;
289
290 /* Flux BCs */
291
292 for (size_t i = 0; i < 3; i++) {
293 af[i] = qimpv[i];
294
295 for (size_t j = 0; j < 3; j++)
296 bf[i][j] = 0;
297 }
298 }
299
300 /*----------------------------------------------------------------------------*/
301 /*!
302 * \brief Set Dirichlet BC for a scalar for a given face.
303 *
304 * \param[out] a explicit BC coefficient for gradients
305 * \param[out] af explicit BC coefficient for diffusive flux
306 * \param[out] b implicit BC coefficient for gradients
307 * \param[out] bf implicit BC coefficient for diffusive flux
308 * \param[in] pimp dirichlet value to impose
309 * \param[in] hint internal exchange coefficient
310 * \param[in] hext external exchange coefficient
311 * (assumed infinite/ignored if < 0)
312 */
313 /*----------------------------------------------------------------------------*/
314
315 inline static void
cs_boundary_conditions_set_dirichlet_scalar(cs_real_t * a,cs_real_t * af,cs_real_t * b,cs_real_t * bf,cs_real_t pimp,cs_real_t hint,cs_real_t hext)316 cs_boundary_conditions_set_dirichlet_scalar(cs_real_t *a,
317 cs_real_t *af,
318 cs_real_t *b,
319 cs_real_t *bf,
320 cs_real_t pimp,
321 cs_real_t hint,
322 cs_real_t hext)
323 {
324 if (hext < 0.) {
325
326 /* Gradient BCs */
327 *a = pimp;
328 *b = 0.;
329
330 /* Flux BCs */
331 *af = -hint*pimp;
332 *bf = hint;
333
334 }
335 else {
336
337 /* Gradient BCs */
338 *a = hext*pimp/(hint + hext);
339 *b = hint /(hint + hext);
340
341 /* Flux BCs */
342 cs_real_t heq = hint*hext/(hint + hext);
343 *af = -heq*pimp;
344 *bf = heq;
345
346 }
347 }
348
349 /*----------------------------------------------------------------------------*/
350 /*!
351 * \brief Set Dirichlet BC for a vector for a given face.
352 *
353 * \param[out] a explicit BC coefficient for gradients
354 * \param[out] af explicit BC coefficient for diffusive flux
355 * \param[out] b implicit BC coefficient for gradients
356 * \param[out] bf implicit BC coefficient for diffusive flux
357 * \param[in] pimpv dirichlet value to impose
358 * \param[in] hint internal exchange coefficient
359 * \param[in] hextv external exchange coefficient
360 * (assumed infinite/ignored if < 0)
361 */
362 /*----------------------------------------------------------------------------*/
363
364 inline static void
cs_boundary_conditions_set_dirichlet_vector(cs_real_3_t a,cs_real_3_t af,cs_real_33_t b,cs_real_33_t bf,cs_real_3_t pimpv,cs_real_t hint,cs_real_3_t hextv)365 cs_boundary_conditions_set_dirichlet_vector(cs_real_3_t a,
366 cs_real_3_t af,
367 cs_real_33_t b,
368 cs_real_33_t bf,
369 cs_real_3_t pimpv,
370 cs_real_t hint,
371 cs_real_3_t hextv)
372 {
373 for (int isou = 0 ; isou < 3; isou++) {
374 if (fabs(hextv[isou]) > 0.5*cs_math_infinite_r) {
375
376 /* Gradient BCs */
377 a[isou] = pimpv[isou];
378 for (int jsou = 0 ; jsou < 3 ; jsou++)
379 b[isou][jsou] = 0.;
380
381 /* Flux BCs */
382 af[isou] = -hint*pimpv[isou];
383 for (int jsou = 0; jsou < 3; jsou++) {
384 if (jsou == isou)
385 bf[isou][jsou] = hint;
386 else
387 bf[isou][jsou] = 0.;
388 }
389 } else {
390
391 cs_real_t heq = hint*hextv[isou]/(hint + hextv[isou]);
392
393 /* Gradient BCs */
394 a[isou] = hextv[isou]*pimpv[isou]/(hint + hextv[isou]);
395 for (int jsou = 0 ; jsou < 3 ; jsou++) {
396 if (jsou == isou)
397 b[isou][jsou] = hint/(hint + hextv[isou]);
398 else
399 b[isou][jsou] = 0.;
400 }
401
402 /* Flux BCs */
403 af[isou] = -heq*pimpv[isou];
404 for (int jsou = 0 ; jsou < 3 ; jsou++) {
405 if (jsou == isou)
406 bf[isou][jsou] = heq;
407 else
408 bf[isou][jsou] = 0.;
409 }
410 }
411 }
412 }
413
414 /*----------------------------------------------------------------------------*/
415 /*!
416 * \brief Set convective oulet boundary condition for a scalar
417 *
418 * Parameters:
419 * \param[out] coefa explicit BC coefficient for gradients
420 * \param[out] cofaf explicit BC coefficient for diffusive flux
421 * \param[out] coefb implicit BC coefficient for gradients
422 * \param[out] cofbf implicit BC coefficient for diffusive flux
423 * \param[in] pimp Flux value to impose
424 * \param[in] cfl Local Courant number used to convect
425 * \param[in] hint Internal exchange coefficient
426 */
427 /*----------------------------------------------------------------------------*/
428
429 void
430 cs_boundary_conditions_set_convective_outlet_scalar(cs_real_t *coefa,
431 cs_real_t *cofaf,
432 cs_real_t *coefb,
433 cs_real_t *cofbf,
434 cs_real_t pimp,
435 cs_real_t cfl,
436 cs_real_t hint);
437
438 /*----------------------------------------------------------------------------*/
439 /*!
440 * \brief Set Dirichlet BC for a vector for a given face with left anisotropic
441 * diffusion.
442 *
443 * \param[out] a explicit BC coefficient for gradients
444 * \param[out] af explicit BC coefficient for diffusive flux
445 * \param[out] b implicit BC coefficient for gradients
446 * \param[out] bf implicit BC coefficient for diffusive flux
447 * \param[in] pimpv dirichlet value to impose
448 * \param[in] hintt internal exchange coefficient
449 * \param[in] hextv external exchange coefficient
450 * (assumed infinite/ignored if < 0)
451 */
452 /*----------------------------------------------------------------------------*/
453
454 inline static void
cs_boundary_conditions_set_dirichlet_vector_aniso(cs_real_3_t a,cs_real_3_t af,cs_real_33_t b,cs_real_33_t bf,cs_real_3_t pimpv,cs_real_6_t hintt,cs_real_3_t hextv)455 cs_boundary_conditions_set_dirichlet_vector_aniso(cs_real_3_t a,
456 cs_real_3_t af,
457 cs_real_33_t b,
458 cs_real_33_t bf,
459 cs_real_3_t pimpv,
460 cs_real_6_t hintt,
461 cs_real_3_t hextv)
462 {
463 for (int isou = 0 ; isou < 3 ; isou++) {
464 if (fabs(hextv[isou]) > 0.5*cs_math_infinite_r) {
465
466 /* Gradient BCs */
467 a[isou] = pimpv[isou];
468 for (int jsou = 0 ; jsou < 3 ; jsou++)
469 b[isou][jsou] = 0.;
470
471 } else {
472
473 cs_exit(1);
474
475 }
476 }
477
478 /* Flux BCs */
479 cs_math_sym_33_3_product(hintt, pimpv, af);
480 for (int isou = 0 ; isou < 3 ; isou++)
481 af[isou] = -af[isou];
482
483 bf[0][0] = hintt[0];
484 bf[1][1] = hintt[1];
485 bf[2][2] = hintt[2];
486 bf[0][1] = hintt[3];
487 bf[1][0] = hintt[3];
488 bf[1][2] = hintt[4];
489 bf[2][1] = hintt[4];
490 bf[0][2] = hintt[5];
491 bf[2][0] = hintt[5];
492 }
493
494 /*----------------------------------------------------------------------------*/
495 /*!
496 * \brief Update per variable boundary condition codes.
497 *
498 * \param[in] nvar number of variables requiring BC's
499 * \param[in] itypfb type of boundary for each face
500 * \param[in, out] icodcl boundary condition codes
501 * \param[in, out] rcodcl boundary condition values
502 */
503 /*----------------------------------------------------------------------------*/
504
505 void
506 cs_boundary_conditions_compute(int nvar,
507 int itypfb[],
508 int icodcl[],
509 double rcodcl[]);
510
511 /*----------------------------------------------------------------------------*/
512 /*!
513 * \brief Automatic adjustments for boundary condition codes.
514 *
515 * Currently handles mapped inlets, after the call to \ref stdtcl.
516 * As portions of stdtcl are migrated to C, they should be called here,
517 * before mapped inlets.
518 *
519 * \param[in] nvar number of variables requiring BC's
520 * \param[in] itypfb type of boundary for each face
521 * \param[in, out] icodcl boundary condition codes
522 * \param[in, out] rcodcl boundary condition values
523 */
524 /*----------------------------------------------------------------------------*/
525
526 void
527 cs_boundary_conditions_complete(int nvar,
528 int itypfb[],
529 int icodcl[],
530 double rcodcl[]);
531
532 /*----------------------------------------------------------------------------*/
533
534 END_C_DECLS
535
536 #endif /* __CS_BOUNDARY_CONDITIONS_H__ */
537