1 /*============================================================================
2 * Porous model management
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 <math.h>
34 #include <float.h>
35 #include <stdio.h>
36 #include <stdlib.h>
37 #include <string.h>
38 #include <assert.h>
39
40 /*----------------------------------------------------------------------------
41 * Local headers
42 *----------------------------------------------------------------------------*/
43
44 #include "bft_mem.h"
45 #include "bft_error.h"
46 #include "bft_printf.h"
47
48 #include "cs_base.h"
49 #include "cs_field.h"
50 #include "cs_log.h"
51 #include "cs_math.h"
52 #include "cs_mesh.h"
53 #include "cs_mesh_quantities.h"
54 #include "cs_preprocess.h"
55
56 /*----------------------------------------------------------------------------
57 * Header for the current file
58 *----------------------------------------------------------------------------*/
59
60 #include "cs_porous_model.h"
61
62 /*----------------------------------------------------------------------------*/
63
64 BEGIN_C_DECLS
65
66 /*----------------------------------------------------------------------------*/
67 /*! \file cs_porous_model.c
68 *
69 * \brief Porous model management
70 */
71 /*----------------------------------------------------------------------------*/
72
73 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
74
75 /*=============================================================================
76 * Local Macro definitions
77 *============================================================================*/
78
79 /*============================================================================
80 * Local type definitions
81 *============================================================================*/
82
83 /*============================================================================
84 * Static global variables
85 *============================================================================*/
86
87 /* Choice of the porous model */
88 int cs_glob_porous_model = 0;
89
90 /*============================================================================
91 * Prototypes for functions intended for use only by Fortran wrappers.
92 * (descriptions follow, with function bodies).
93 *============================================================================*/
94
95 void
96 cs_f_porous_model_get_pointers(int **iporos);
97
98 int
99 cs_f_porous_model_cell_is_active(cs_lnum_t cell_id);
100
101 /*=============================================================================
102 * Private function definitions
103 *============================================================================*/
104
105 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
106
107 /*============================================================================
108 * Fortran wrapper function definitions
109 *============================================================================*/
110
111 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
112
113 /*----------------------------------------------------------------------------
114 * Get pointers to global variables.
115 *
116 * This function is intended for use by Fortran wrappers, and
117 * enables mapping to Fortran global pointers.
118 *
119 * parameters:
120 * iporos --> pointer to cs_glob_porous_model
121 *----------------------------------------------------------------------------*/
122
123 void
cs_f_porous_model_get_pointers(int ** iporos)124 cs_f_porous_model_get_pointers(int **iporos)
125 {
126 *iporos = &(cs_glob_porous_model);
127 }
128
129 /*----------------------------------------------------------------------------*/
130 /*!
131 * \brief Return 0 if cell is disabled, 1 otherwise
132 *
133 * \param[in] cell_id
134 *
135 * \return 0 if cell is disabled
136 */
137 /*----------------------------------------------------------------------------*/
138
139 int
cs_f_porous_model_cell_is_active(cs_lnum_t cell_id)140 cs_f_porous_model_cell_is_active(cs_lnum_t cell_id)
141 {
142 cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
143
144 return (1 - (mq->has_disable_flag
145 * mq->c_disable_flag[mq->has_disable_flag * cell_id]));
146 }
147
148 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
149
150 /*=============================================================================
151 * Public function definitions
152 *============================================================================*/
153
154 /*----------------------------------------------------------------------------
155 * Set porous model option.
156 *
157 * parameters:
158 * porous_model <-- porous model option (> 0 for porosity)
159 *----------------------------------------------------------------------------*/
160
161 void
cs_porous_model_set_model(int porous_model)162 cs_porous_model_set_model(int porous_model)
163 {
164 cs_glob_porous_model = porous_model;
165 }
166
167 /*----------------------------------------------------------------------------*/
168 /*!
169 * \brief Initialize disable_flag
170 */
171 /*----------------------------------------------------------------------------*/
172
173 void
cs_porous_model_init_disable_flag(void)174 cs_porous_model_init_disable_flag(void)
175 {
176 cs_mesh_t *m =cs_glob_mesh;
177 cs_mesh_quantities_t *mq =cs_glob_mesh_quantities;
178
179 static cs_lnum_t n_cells_prev = 0;
180
181 if (cs_glob_porous_model > 0) {
182 cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
183 if (mq->c_disable_flag == NULL || m->n_cells != n_cells_prev) {
184 /* Currently should handle turbomachinery (ghost cells only can change),
185 not mesh refinement */
186 assert(m->n_cells == n_cells_prev || n_cells_prev == 0);
187
188 BFT_REALLOC(mq->c_disable_flag, n_cells_ext, int);
189 for (cs_lnum_t cell_id = 0; cell_id < n_cells_ext; cell_id++)
190 mq->c_disable_flag[cell_id] = 0;
191
192 n_cells_prev = m->n_cells;
193 }
194 else {
195 if (mq->has_disable_flag != 0)
196 BFT_REALLOC(mq->c_disable_flag, n_cells_ext, int);
197 if (m->halo != NULL)
198 cs_halo_sync_untyped(m->halo, CS_HALO_STANDARD, sizeof(int),
199 mq->c_disable_flag);
200 }
201 }
202 else {
203 if (mq->c_disable_flag == NULL)
204 BFT_MALLOC(mq->c_disable_flag, 1, int);
205 mq->c_disable_flag[0] = 0;
206 }
207
208 /* Update Fortran pointer quantities */
209 cs_preprocess_mesh_update_fortran();
210 }
211
212 /*----------------------------------------------------------------------------*/
213 /*!
214 * \brief Set (unset) has_disable_flag
215 *
216 * \param[in] flag 1: on, 0: off
217 */
218 /*----------------------------------------------------------------------------*/
219
220 void
cs_porous_model_set_has_disable_flag(int flag)221 cs_porous_model_set_has_disable_flag(int flag)
222 {
223 cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
224
225 mq->has_disable_flag = flag;
226
227 /* if off, fluid surfaces point toward cell surfaces */
228 /* Porous models */
229 if (cs_glob_porous_model > 0) {
230 if (flag == 0) {
231 /* Set pointers of face quantities to the standard ones */
232 if (cs_glob_porous_model == 3) {
233 mq->i_f_face_normal = mq->i_face_normal;
234 mq->b_f_face_normal = mq->b_face_normal;
235 mq->i_f_face_surf = mq->i_face_surf;
236 mq->b_f_face_surf = mq->b_face_surf;
237 mq->i_f_face_factor = NULL;
238 mq->b_f_face_factor = NULL;
239 }
240 mq->cell_f_vol = mq->cell_vol;
241 }
242 else {
243 /* Use fluid surfaces and volumes */
244 if (cs_glob_porous_model == 3) {
245 mq->i_f_face_normal = cs_field_by_name("i_f_face_normal")->val;
246 mq->b_f_face_normal = cs_field_by_name("b_f_face_normal")->val;
247 mq->i_f_face_surf = cs_field_by_name("i_f_face_surf")->val;
248 mq->b_f_face_surf = cs_field_by_name("b_f_face_surf")->val;
249 mq->i_f_face_factor
250 = (cs_real_2_t *)cs_field_by_name("i_f_face_factor")->val;
251 mq->b_f_face_factor = cs_field_by_name("b_f_face_factor")->val;
252 }
253 mq->cell_f_vol = cs_field_by_name("cell_f_vol")->val;
254 }
255 }
256
257 /* Update Fortran pointer quantities */
258 cs_preprocess_mesh_update_fortran();
259 }
260
261 /*----------------------------------------------------------------------------*/
262 /*!
263 * \brief Init fluid quantities
264 */
265 /*----------------------------------------------------------------------------*/
266
267 void
cs_porous_model_init_fluid_quantities(void)268 cs_porous_model_init_fluid_quantities(void)
269 {
270 if (cs_glob_porous_model == 3) {
271 cs_mesh_init_fluid_sections(cs_glob_mesh, cs_glob_mesh_quantities);
272 }
273 }
274
275 /*----------------------------------------------------------------------------*/
276 /*!
277 * \brief Automatic computation of the face porosity and factors.
278 *
279 * This is useful for the integral porous model.
280 */
281 /*----------------------------------------------------------------------------*/
282
283 void
cs_porous_model_auto_face_porosity(void)284 cs_porous_model_auto_face_porosity(void)
285 {
286 if (cs_glob_porous_model < 3)
287 return;
288
289 cs_mesh_t *m = cs_glob_mesh;
290 cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
291
292 /* Get the cell porosity field value */
293 cs_real_t *cpro_porosi = cs_field_by_name("porosity")->val;
294
295 if (m->halo != NULL)
296 cs_halo_sync_var(m->halo, CS_HALO_STANDARD, cpro_porosi);
297
298 /* Set interior face values */
299
300 {
301 const cs_lnum_t n_i_faces = m->n_i_faces;
302 const cs_lnum_2_t *i_face_cells
303 = (const cs_lnum_2_t *)m->i_face_cells;
304
305 const cs_real_3_t *restrict i_face_normal
306 = (const cs_real_3_t *restrict)mq->i_face_normal;
307 cs_real_3_t *restrict i_f_face_normal
308 = (cs_real_3_t *restrict)mq->i_f_face_normal;
309
310 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
311
312 cs_lnum_t c_id0 = i_face_cells[face_id][0];
313 cs_lnum_t c_id1 = i_face_cells[face_id][1];
314
315 cs_real_t face_porosity = CS_MIN(cpro_porosi[c_id0], cpro_porosi[c_id1]);
316
317 for (cs_lnum_t i = 0; i < 3; i++)
318 i_f_face_normal[face_id][i] = face_porosity * i_face_normal[face_id][i];
319
320 mq->i_f_face_surf[face_id] = cs_math_3_norm(i_f_face_normal[face_id]);
321
322 if (mq->i_f_face_factor != NULL) {
323 if (face_porosity > cs_math_epzero) {
324 mq->i_f_face_factor[face_id][0] = cpro_porosi[c_id0] / face_porosity;
325 mq->i_f_face_factor[face_id][1] = cpro_porosi[c_id1] / face_porosity;
326 }
327 else {
328 mq->i_f_face_factor[face_id][0] = 1.;
329 mq->i_f_face_factor[face_id][1] = 1;
330 }
331 }
332
333 }
334 }
335
336 /* Set boundary face values */
337
338 {
339 const cs_lnum_t n_b_faces = m->n_b_faces;
340 const cs_lnum_t *b_face_cells
341 = (const cs_lnum_t *)m->b_face_cells;
342
343 const cs_real_3_t *restrict b_face_normal
344 = (const cs_real_3_t *restrict)mq->b_face_normal;
345 cs_real_3_t *restrict b_f_face_normal
346 = (cs_real_3_t *restrict)mq->b_f_face_normal;
347
348 for (cs_lnum_t face_id = 0; face_id < n_b_faces; face_id++) {
349
350 cs_lnum_t c_id = b_face_cells[face_id];
351
352 cs_real_t face_porosity = cpro_porosi[c_id];
353
354 for (cs_lnum_t i = 0; i < 3; i++)
355 b_f_face_normal[face_id][i] = face_porosity * b_face_normal[face_id][i];
356
357 mq->b_f_face_surf[face_id] = cs_math_3_norm(b_f_face_normal[face_id]);
358
359 if (mq->b_f_face_factor != NULL) {
360 if (face_porosity > cs_math_epzero) {
361 mq->b_f_face_factor[face_id] = cpro_porosi[c_id] / face_porosity;
362 }
363 else {
364 mq->b_f_face_factor[face_id] = 1.;
365 }
366 }
367
368 }
369 }
370 }
371
372 /*----------------------------------------------------------------------------*/
373
374 END_C_DECLS
375