1 /*============================================================================
2 * User definitions of porous media.
3 *============================================================================*/
4
5 /* VERS */
6
7 /*
8 This file is part of Code_Saturne, a general-purpose CFD tool.
9
10 Copyright (C) 1998-2021 EDF S.A.
11
12 This program is free software; you can redistribute it and/or modify it under
13 the terms of the GNU General Public License as published by the Free Software
14 Foundation; either version 2 of the License, or (at your option) any later
15 version.
16
17 This program is distributed in the hope that it will be useful, but WITHOUT
18 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
20 details.
21
22 You should have received a copy of the GNU General Public License along with
23 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
24 Street, Fifth Floor, Boston, MA 02110-1301, USA.
25 */
26
27 /*----------------------------------------------------------------------------*/
28
29 #include "cs_defs.h"
30
31 /*----------------------------------------------------------------------------
32 * Standard C library headers
33 *----------------------------------------------------------------------------*/
34
35 #include <assert.h>
36 #include <math.h>
37 #include <string.h>
38
39 #if defined(HAVE_MPI)
40 #include <mpi.h>
41 #endif
42
43 /*----------------------------------------------------------------------------
44 * Local headers
45 *----------------------------------------------------------------------------*/
46
47 #include "cs_headers.h"
48
49 /*----------------------------------------------------------------------------*/
50
51 BEGIN_C_DECLS
52
53 /*----------------------------------------------------------------------------*/
54 /*!
55 * \file cs_user_porosity.c
56 *
57 * \brief User definitions of porous media.
58 *
59 * See \ref cs_porosity for examples.
60 */
61 /*----------------------------------------------------------------------------*/
62
63 /*============================================================================
64 * User function definitions
65 *============================================================================*/
66
67 /*----------------------------------------------------------------------------*/
68 /*!
69 * \brief Compute the porosity (volume factor \f$ \epsilon \f$
70 * when the porosity model is activated.
71 * (\ref cs_glob_porous_model > 0).
72 *
73 * This function is called at the beginning of the simulation only.
74 *
75 * \param[in, out] domain pointer to a cs_domain_t structure
76 */
77 /*----------------------------------------------------------------------------*/
78
79 void
cs_user_porosity(cs_domain_t * domain)80 cs_user_porosity(cs_domain_t *domain)
81 {
82 /*!< [init_poro_mq] */
83 cs_mesh_t *m = domain->mesh;
84 cs_mesh_quantities_t *mq = domain->mesh_quantities;
85
86 const cs_lnum_2_t *i_face_cells
87 = (const cs_lnum_2_t *)m->i_face_cells;
88
89 const cs_real_3_t *restrict i_face_cog
90 = (const cs_real_3_t *restrict)mq->i_face_cog;
91 const cs_real_3_t *restrict b_face_cog
92 = (const cs_real_3_t *restrict)mq->b_face_cog;
93 const cs_real_3_t *restrict cell_cen
94 = (const cs_real_3_t *restrict)mq->cell_cen;
95 const cs_real_3_t *restrict i_face_normal
96 = (const cs_real_3_t *restrict)mq->i_face_normal;
97 cs_real_3_t *restrict i_f_face_normal
98 = (cs_real_3_t *restrict)mq->i_f_face_normal;
99 const cs_real_3_t *restrict b_face_normal
100 = (const cs_real_3_t *restrict)mq->b_face_normal;
101 cs_real_3_t *restrict b_f_face_normal
102 = (cs_real_3_t *restrict)mq->b_f_face_normal;
103
104 const cs_real_t *i_f_face_surf = mq->i_f_face_surf;
105 const cs_real_t *i_face_surf = mq->i_face_surf;
106 /*!< [init_poro_mq] */
107
108 /* Get the cell porosity field value */
109 /*!< [init_poro_pro] */
110 cs_real_t *cpro_porosi = cs_field_by_name("porosity")->val;
111 /*!< [init_poro_pro] */
112
113 /* First, set cell porosity value; the fluid cell volume will be
114 * automatically deduced */
115
116 /*!< [set_poro_cells_1] */
117 for (cs_lnum_t cell_id = 0; cell_id < m->n_cells; cell_id++) {
118
119 cs_real_t x = cell_cen[cell_id][0];
120
121 if (x < 20.)
122 cpro_porosi[cell_id] = 1.;
123 else
124 cpro_porosi[cell_id] = 0.5;
125 }
126
127 /* synchronize for use in fluid face factor calculation below */
128 cs_halo_type_t halo_type = CS_HALO_STANDARD;
129 cs_field_synchronize(cs_field_by_name("porosity"),
130 halo_type);
131
132 /*!< [set_poro_cells_1] */
133
134 /* Set interior face values (for cs_glob_porous_model == 3 only) */
135
136 /*!< [set_poro_i_faces_1] */
137 for (cs_lnum_t face_id = 0; face_id < m->n_i_faces; face_id++) {
138
139 cs_real_t x = i_face_cog[face_id][0];
140
141 cs_real_t face_porosity = 1.;
142 if (x > 19.9999)
143 face_porosity = 0.5;
144
145 for (int i = 0; i < 3; i++)
146 i_f_face_normal[face_id][i] = face_porosity * i_face_normal[face_id][i];
147
148 mq->i_f_face_surf[face_id] = cs_math_3_norm(i_f_face_normal[face_id]);
149
150 }
151 /*!< [set_poro_i_faces_1] */
152
153 /* Set boundary face values (for cs_glob_porous_model == 3 only) */
154
155 /*!< [set_poro_b_faces_1] */
156 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
157
158 cs_real_t x = b_face_cog[face_id][0];
159
160 cs_real_t face_porosity = 1.;
161 if (x > 19.9999)
162 face_porosity = 0.5;
163
164 for (int i = 0; i < 3; i++)
165 b_f_face_normal[face_id][i] = face_porosity * b_face_normal[face_id][i];
166
167 mq->b_f_face_surf[face_id] = cs_math_3_norm(b_f_face_normal[face_id]);
168
169 }
170 /*!< [set_poro_b_faces_1] */
171
172 /* Four set face factor */
173 if (mq->i_f_face_factor != NULL) {
174 for (cs_lnum_t face_id = 0; face_id < m->n_i_faces; face_id++) {
175 cs_lnum_t ii = i_face_cells[face_id][0];
176 cs_lnum_t jj = i_face_cells[face_id][1];
177
178 cs_real_t face_porosity =
179 CS_MAX(
180 i_f_face_surf[face_id] / i_face_surf[face_id],
181 cs_math_epzero);
182
183 mq->i_f_face_factor[face_id][0] = cpro_porosi[ii] / face_porosity;
184 mq->i_f_face_factor[face_id][1] = cpro_porosi[jj] / face_porosity;
185 }
186 }
187
188 }
189
190 /*----------------------------------------------------------------------------*/
191
192 END_C_DECLS
193