1 /*============================================================================
2 * Filters for dynamic models.
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 <errno.h>
35 #include <stdio.h>
36 #include <stdarg.h>
37 #include <string.h>
38 #include <math.h>
39 #include <float.h>
40
41 #if defined(HAVE_MPI)
42 #include <mpi.h>
43 #endif
44
45 /*----------------------------------------------------------------------------
46 * Local headers
47 *----------------------------------------------------------------------------*/
48
49 #include "bft_error.h"
50 #include "bft_mem.h"
51 #include "bft_printf.h"
52
53 #include "cs_cell_to_vertex.h"
54 #include "cs_ext_neighborhood.h"
55 #include "cs_mesh.h"
56 #include "cs_mesh_adjacencies.h"
57 #include "cs_mesh_quantities.h"
58
59 /*----------------------------------------------------------------------------
60 * Header for the current file
61 *----------------------------------------------------------------------------*/
62
63 #include "cs_les_filter.h"
64 #include "cs_halo.h"
65 #include "cs_halo_perio.h"
66
67 /*----------------------------------------------------------------------------*/
68
69 BEGIN_C_DECLS
70
71 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
72
73 /*=============================================================================
74 * Local Macro definitions
75 *============================================================================*/
76
77 /*============================================================================
78 * Type definition
79 *============================================================================*/
80
81 /*============================================================================
82 * Static global variables
83 *============================================================================*/
84
85 /*============================================================================
86 * Private function definitions
87 *============================================================================*/
88
89 /*----------------------------------------------------------------------------*/
90 /*!
91 * \brief Compute filters for dynamic models.
92 *
93 * This function deals with the standard or extended neighborhood.
94 *
95 * \param[in] stride stride of array to filter
96 * \param[in] val array of values to filter
97 * \param[out] f_val array of filtered values
98 */
99 /*----------------------------------------------------------------------------*/
100
101 static void
_les_filter_ext_neighborhood(int stride,cs_real_t val[],cs_real_t f_val[])102 _les_filter_ext_neighborhood(int stride,
103 cs_real_t val[],
104 cs_real_t f_val[])
105 {
106 cs_real_t *w1 = NULL, *w2 = NULL;
107
108 const cs_mesh_t *mesh = cs_glob_mesh;
109 const cs_lnum_t _stride = stride;
110 const cs_lnum_t n_cells = mesh->n_cells;
111 const cs_lnum_t n_cells_ext = mesh->n_cells_with_ghosts;
112 const cs_lnum_t n_elts_l = n_cells * _stride;
113 const cs_lnum_t n_elts = n_cells_ext * _stride;
114 const int n_i_groups = mesh->i_face_numbering->n_groups;
115 const int n_i_threads = mesh->i_face_numbering->n_threads;
116 const cs_lnum_t *restrict i_group_index = mesh->i_face_numbering->group_index;
117 const cs_lnum_t *cell_cells_idx = mesh->cell_cells_idx;
118 const cs_lnum_t *cell_cells_lst = mesh->cell_cells_lst;
119 const cs_real_t *cell_vol = cs_glob_mesh_quantities->cell_vol;
120
121 assert(cell_cells_idx != NULL);
122
123 /* Allocate and initialize working buffers */
124
125 BFT_MALLOC(w1, n_elts, cs_real_t);
126 BFT_MALLOC(w2, n_elts, cs_real_t);
127
128 /* Case for scalar variable */
129 /*--------------------------*/
130
131 if (stride == 1) {
132
133 /* Synchronize variable */
134
135 if (mesh->halo != NULL)
136 cs_halo_sync_var(mesh->halo, CS_HALO_EXTENDED, val);
137
138 /* Define filtered variable array */
139
140 # pragma omp parallel for
141 for (cs_lnum_t i = 0; i < n_cells; i++) {
142
143 w1[i] = val[i] * cell_vol[i];
144 w2[i] = cell_vol[i];
145
146 /* Loop on connected cells (without cells sharing a face) */
147
148 for (cs_lnum_t j = cell_cells_idx[i]; j < cell_cells_idx[i+1]; j++) {
149 cs_lnum_t k = cell_cells_lst[j];
150 w1[i] += val[k] * cell_vol[k];
151 w2[i] += cell_vol[k];
152 }
153
154 } /* End of loop on cells */
155
156 # pragma omp parallel for if (mesh->n_ghost_cells > CS_THR_MIN)
157 for (cs_lnum_t i = n_cells; i < n_cells_ext; i++) {
158 w1[i] = 0;
159 w2[i] = 0;
160 }
161
162 for (int g_id = 0; g_id < n_i_groups; g_id++) {
163
164 # pragma omp parallel for
165 for (int t_id = 0; t_id < n_i_threads; t_id++) {
166
167 for (cs_lnum_t face_id = i_group_index[(t_id*n_i_groups + g_id)*2];
168 face_id < i_group_index[(t_id*n_i_groups + g_id)*2 + 1];
169 face_id++) {
170
171 cs_lnum_t i = mesh->i_face_cells[face_id][0];
172 cs_lnum_t j = mesh->i_face_cells[face_id][1];
173
174 w1[i] += val[j] * cell_vol[j];
175 w2[i] += cell_vol[j];
176
177 w1[j] += val[i] * cell_vol[i];
178 w2[j] += cell_vol[i];
179 }
180
181 }
182
183 }
184
185 # pragma omp parallel for
186 for (cs_lnum_t i = 0; i < n_cells; i++)
187 f_val[i] = w1[i]/w2[i];
188
189 /* Synchronize variable */
190
191 if (mesh->halo != NULL)
192 cs_halo_sync_var(mesh->halo, CS_HALO_STANDARD, f_val);
193 }
194
195 /* Case for vector or tensor variable */
196 /*------------------------------------*/
197
198 else {
199
200 /* Synchronize variable */
201
202 if (mesh->halo != NULL)
203 cs_halo_sync_var_strided(mesh->halo, CS_HALO_EXTENDED, val, stride);
204
205 /* Define filtered variable array */
206
207 # pragma omp parallel for
208 for (cs_lnum_t i = 0; i < n_cells; i++) {
209
210 for (cs_lnum_t c_id = 0; c_id < _stride; c_id++) {
211 const cs_lnum_t ic = i *_stride + c_id;
212 w1[ic] = val[ic] * cell_vol[i];
213 w2[ic] = cell_vol[i];
214 }
215
216 /* Loop on connected cells (without cells sharing a face) */
217
218 for (cs_lnum_t j = cell_cells_idx[i]; j < cell_cells_idx[i+1]; j++) {
219 cs_lnum_t k = cell_cells_lst[j];
220 for (cs_lnum_t c_id = 0; c_id < _stride; c_id++) {
221 const cs_lnum_t ic = i *_stride + c_id;
222 const cs_lnum_t kc = k *_stride + c_id;
223 w1[ic] += val[kc] * cell_vol[k];
224 w2[ic] += cell_vol[k];
225 }
226 }
227
228 } /* End of loop on cells */
229
230 # pragma omp parallel for if (mesh->n_ghost_cells > CS_THR_MIN)
231 for (cs_lnum_t i = n_cells; i < n_cells_ext; i++) {
232 for (cs_lnum_t c_id = 0; c_id < _stride; c_id++) {
233 const cs_lnum_t ic = i *_stride + c_id;
234 w1[ic] = 0;
235 w2[ic] = 0;
236 }
237 }
238
239 for (int g_id = 0; g_id < n_i_groups; g_id++) {
240
241 # pragma omp parallel for
242 for (int t_id = 0; t_id < n_i_threads; t_id++) {
243
244 for (cs_lnum_t face_id = i_group_index[(t_id*n_i_groups + g_id)*2];
245 face_id < i_group_index[(t_id*n_i_groups + g_id)*2 + 1];
246 face_id++) {
247
248 cs_lnum_t i = mesh->i_face_cells[face_id][0];
249 cs_lnum_t j = mesh->i_face_cells[face_id][1];
250
251 for (cs_lnum_t c_id = 0; c_id < _stride; c_id++) {
252 const cs_lnum_t ic = i *_stride + c_id;
253 const cs_lnum_t jc = j *_stride + c_id;
254
255 w1[ic] += val[jc] * cell_vol[j];
256 w2[ic] += cell_vol[j];
257
258 w1[jc] += val[ic] * cell_vol[i];
259 w2[jc] += cell_vol[i];
260 }
261
262 }
263
264 }
265
266 }
267
268 # pragma omp parallel for
269 for (cs_lnum_t i = 0; i < n_elts_l; i++)
270 f_val[i] = w1[i]/w2[i];
271
272 /* Synchronize variable */
273
274 if (mesh->halo != NULL)
275 cs_halo_sync_var_strided(mesh->halo, CS_HALO_EXTENDED, f_val, stride);
276 }
277
278 BFT_FREE(w2);
279 BFT_FREE(w1);
280 }
281
282 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
283
284 /*=============================================================================
285 * Public function definitions
286 *============================================================================*/
287
288 /*----------------------------------------------------------------------------*/
289 /*!
290 * \brief Compute filters for dynamic models.
291 *
292 * This function deals with the standard or extended neighborhood.
293 *
294 * \param[in] stride stride of array to filter
295 * \param[in] val array of values to filter
296 * \param[out] f_val array of filtered values
297 */
298 /*----------------------------------------------------------------------------*/
299
300 void
cs_les_filter(int stride,cs_real_t val[],cs_real_t f_val[])301 cs_les_filter(int stride,
302 cs_real_t val[],
303 cs_real_t f_val[])
304 {
305 if (cs_ext_neighborhood_get_type() == CS_EXT_NEIGHBORHOOD_COMPLETE) {
306 _les_filter_ext_neighborhood(stride, val, f_val);
307 return;
308 }
309
310 cs_real_t *v_val = NULL, *v_weight = NULL;
311
312 const cs_mesh_t *mesh = cs_glob_mesh;
313 const cs_lnum_t _stride = stride;
314 const cs_lnum_t n_cells = mesh->n_cells;
315 const cs_real_t *cell_vol = cs_glob_mesh_quantities->cell_vol;
316
317 /* Allocate and initialize working buffer */
318
319 BFT_MALLOC(v_val, mesh->n_vertices*_stride, cs_real_t);
320 BFT_MALLOC(v_weight, mesh->n_vertices, cs_real_t);
321
322 /* Define filtered variable array */
323
324 cs_cell_to_vertex(CS_CELL_TO_VERTEX_LR,
325 0,
326 _stride,
327 true, /* ignore periodicity of rotation */
328 cell_vol,
329 val,
330 NULL,
331 v_val);
332
333 cs_cell_to_vertex(CS_CELL_TO_VERTEX_LR,
334 0,
335 1,
336 true, /* ignore periodicity of rotation */
337 NULL,
338 cell_vol,
339 NULL,
340 v_weight);
341
342 /* Build cell average */
343
344 const cs_adjacency_t *c2v = cs_mesh_adjacencies_cell_vertices();
345 const cs_lnum_t *c2v_idx = c2v->idx;
346 const cs_lnum_t *c2v_ids = c2v->ids;
347
348 if (_stride == 1) {
349
350 # pragma omp parallel for if (n_cells > CS_THR_MIN)
351 for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
352 cs_lnum_t s_id = c2v_idx[c_id];
353 cs_lnum_t e_id = c2v_idx[c_id+1];
354 cs_real_t _f_val = 0, _f_weight = 0;
355 for (cs_lnum_t j = s_id; j < e_id; j++) {
356 cs_lnum_t v_id = c2v_ids[j];
357 _f_val += v_val[v_id] * v_weight[v_id];
358 _f_weight += v_weight[v_id];
359 }
360 f_val[c_id] = _f_val / _f_weight;
361 }
362
363 }
364 else {
365
366 assert(_stride <= 9);
367
368 # pragma omp parallel for if (n_cells > CS_THR_MIN)
369 for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
370 cs_lnum_t s_id = c2v_idx[c_id];
371 cs_lnum_t e_id = c2v_idx[c_id+1];
372 cs_real_t _f_val[9], _f_weight = 0;
373 for (cs_lnum_t k = 0; k < _stride; k++)
374 _f_val[k] = 0;
375 for (cs_lnum_t j = s_id; j < e_id; j++) {
376 cs_lnum_t v_id = c2v_ids[j];
377 for (cs_lnum_t k = 0; k < _stride; k++) {
378 _f_val[k] += v_val[v_id*_stride + k] * v_weight[v_id];
379 }
380 _f_weight += v_weight[v_id];
381 }
382 for (cs_lnum_t k = 0; k < _stride; k++) {
383 f_val[c_id*_stride + k] = _f_val[k] / _f_weight;
384 }
385
386 }
387
388 }
389
390 BFT_FREE(v_weight);
391 BFT_FREE(v_val);
392
393 /* Synchronize variable */
394
395 if (mesh->halo != NULL)
396 cs_halo_sync_var_strided(mesh->halo, CS_HALO_STANDARD, f_val, _stride);
397 }
398
399 /*----------------------------------------------------------------------------*/
400
401 END_C_DECLS
402