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