1 /*============================================================================
2 * Regulation on bad cells.
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 <stdio.h>
34 #include <stdlib.h>
35 #include <string.h>
36 #include <math.h>
37
38 /*----------------------------------------------------------------------------
39 * BFT library headers
40 *----------------------------------------------------------------------------*/
41
42 #include "bft_mem.h"
43 #include "bft_printf.h"
44
45 /*----------------------------------------------------------------------------
46 * Code_Saturne library headers
47 *----------------------------------------------------------------------------*/
48
49 #include "cs_blas.h"
50 #include "cs_boundary_conditions.h"
51 #include "cs_halo.h"
52 #include "cs_halo_perio.h"
53 #include "cs_mesh.h"
54 #include "cs_mesh_bad_cells.h"
55 #include "cs_mesh_quantities.h"
56 #include "cs_parall.h"
57 #include "cs_parameters.h"
58 #include "cs_sles_default.h"
59
60 /*----------------------------------------------------------------------------
61 * Header for the current file
62 *----------------------------------------------------------------------------*/
63
64 #include "cs_bad_cells_regularisation.h"
65
66 /*----------------------------------------------------------------------------*/
67
68 #ifndef DOXYGEN_SHOULD_SKIP_THIS
69
70 BEGIN_C_DECLS
71
72 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
73
74 /*----------------------------------------------------------------------------*/
75 /*!
76 * \brief Regularisation on bad cells for scalars
77 *
78 */
79 /*----------------------------------------------------------------------------*/
80
81 void
cs_bad_cells_regularisation_scalar(cs_real_t * var)82 cs_bad_cells_regularisation_scalar(cs_real_t *var)
83 {
84 const cs_mesh_t *mesh = cs_glob_mesh;
85 cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
86
87 if (!(cs_glob_mesh_quantities_flag & CS_BAD_CELLS_REGULARISATION))
88 return;
89
90 cs_lnum_t n_cells_ext = mesh->n_cells_with_ghosts;
91 cs_lnum_t n_cells = mesh->n_cells;
92 cs_lnum_t n_i_faces = mesh->n_i_faces;
93 const cs_lnum_2_t *i_face_cells = (const cs_lnum_2_t *)mesh->i_face_cells;
94
95 const cs_real_t *surfn = mq->i_face_surf;
96 double *dist = mq->i_dist;
97 double *volume = mq->cell_vol;
98
99 cs_real_t *xam, *dam, *rhs;
100
101 double varmin = 1.e20;
102 double varmax =-1.e20;
103
104 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++)
105 if (!(mq->bad_cell_flag[cell_id] & CS_BAD_CELL_TO_REGULARIZE)) {
106 varmin = CS_MIN(varmin, var[cell_id]);
107 varmax = CS_MAX(varmax, var[cell_id]);
108 }
109
110 cs_parall_min(1, CS_DOUBLE, &varmin);
111 cs_parall_max(1, CS_DOUBLE, &varmax);
112
113 BFT_MALLOC(xam, n_i_faces, cs_real_t);
114 BFT_MALLOC(dam, n_cells_ext, cs_real_t);
115 BFT_MALLOC(rhs, n_cells_ext, cs_real_t);
116
117 for (cs_lnum_t cell_id = 0; cell_id < n_cells_ext; cell_id++) {
118 dam[cell_id] = 0.;
119 rhs[cell_id] = 0.;
120 }
121
122 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
123 cs_lnum_t cell_id1 = i_face_cells[face_id][0];
124 cs_lnum_t cell_id2 = i_face_cells[face_id][1];
125 xam[face_id] = 0.;
126
127 double surf = surfn[face_id];
128 double vol = 0.5 * (volume[cell_id1] + volume[cell_id2]);
129 surf = CS_MAX(surf, 0.1*vol/dist[face_id]);
130 double ssd = surf / dist[face_id];
131
132 dam[cell_id1] += ssd;
133 dam[cell_id2] += ssd;
134
135 if ( mq->bad_cell_flag[cell_id1] & CS_BAD_CELL_TO_REGULARIZE
136 && mq->bad_cell_flag[cell_id2] & CS_BAD_CELL_TO_REGULARIZE) {
137 xam[face_id] = -ssd;
138 }
139 else if (mq->bad_cell_flag[cell_id1] & CS_BAD_CELL_TO_REGULARIZE) {
140 rhs[cell_id1] += ssd * var[cell_id2];
141 rhs[cell_id2] += ssd * var[cell_id2];
142 }
143 else if (mq->bad_cell_flag[cell_id2] & CS_BAD_CELL_TO_REGULARIZE) {
144 rhs[cell_id2] += ssd * var[cell_id1];
145 rhs[cell_id1] += ssd * var[cell_id1];
146 }
147 else {
148 rhs[cell_id1] += ssd * var[cell_id1];
149 rhs[cell_id2] += ssd * var[cell_id2];
150 }
151 }
152
153 cs_real_t rnorm = sqrt(cs_gdot(n_cells, rhs, rhs));
154
155 /* Solver residual */
156 double ressol = 0.;
157 int niterf = 0;
158
159 cs_real_t epsilp = 1.e-12;
160
161 /* Matrix block size */
162 cs_lnum_t db_size[4] = {1, 1, 1, 1};
163
164 cs_sles_solve_native(-1, /* f_id */
165 "potential_regularisation_scalar",
166 true, /* symmetric */
167 db_size,
168 NULL, /* eb_size */
169 (cs_real_t *)dam,
170 xam,
171 epsilp,
172 rnorm,
173 &niterf,
174 &ressol,
175 (cs_real_t *)rhs,
176 (cs_real_t *)var);
177
178 bft_printf("Solving %s: N iter: %d, Res: %12.5e, Norm: %12.5e\n",
179 "potential_regularisation_scalar", niterf, ressol, rnorm);
180 //FIXME use less!
181 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
182 var[cell_id] = CS_MIN(var[cell_id], varmax);
183 var[cell_id] = CS_MAX(var[cell_id], varmin);
184 }
185
186 if (mesh->halo != NULL)
187 cs_halo_sync_var(mesh->halo, CS_HALO_STANDARD, var);
188
189 /* Free solver setup */
190 cs_sles_free_native(-1, /* f_id*/
191 "potential_regularisation_scalar");
192
193 /* Free memory */
194 BFT_FREE(xam);
195 BFT_FREE(dam);
196 BFT_FREE(rhs);
197
198 return;
199 }
200
201 /*----------------------------------------------------------------------------*/
202 /*!
203 * \brief Regularisation on bad cells for vectors
204 *
205 */
206 /*----------------------------------------------------------------------------*/
207
208 void
cs_bad_cells_regularisation_vector(cs_real_3_t * var,int boundary_projection)209 cs_bad_cells_regularisation_vector(cs_real_3_t *var,
210 int boundary_projection)
211 {
212 const cs_mesh_t *mesh = cs_glob_mesh;
213 cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
214
215 if (!(cs_glob_mesh_quantities_flag & CS_BAD_CELLS_REGULARISATION))
216 return;
217
218 cs_lnum_t n_cells_ext = mesh->n_cells_with_ghosts;
219 cs_lnum_t n_cells = mesh->n_cells;
220 cs_lnum_t n_i_faces = mesh->n_i_faces;
221 cs_lnum_t n_b_faces = mesh->n_b_faces;
222 const cs_lnum_2_t *i_face_cells = (const cs_lnum_2_t *)mesh->i_face_cells;
223 const cs_lnum_t *b_face_cells = mesh->b_face_cells;
224
225 const cs_real_t *surfn = mq->i_face_surf;
226 const cs_real_t *surfbn = mq->b_face_surf;
227 double *dist = mq->i_dist;
228 double *distbr = mq->b_dist;
229 double *volume = mq->cell_vol;
230
231 const cs_real_3_t *surfbo = (const cs_real_3_t *) mq->b_face_normal;
232
233 cs_real_33_t *dam;
234 cs_real_3_t *rhs;
235 cs_real_t *xam;
236 #if 1
237 double varmin[3] = {1.e20, 1.e20, 1.e20};
238 double varmax[3] = {-1.e20, -1.e20,-1.e20};
239
240 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++)
241 if (!(mq->bad_cell_flag[cell_id] & CS_BAD_CELL_TO_REGULARIZE)) {
242 for (int i = 0; i < 3; i++) {
243 varmin[i] = CS_MIN(varmin[i], var[cell_id][i]);
244 varmax[i] = CS_MAX(varmax[i], var[cell_id][i]);
245 }
246 }
247
248 for (int i = 0; i < 3; i++) {
249 cs_parall_min(1, CS_DOUBLE, &varmin[i]);
250 cs_parall_max(1, CS_DOUBLE, &varmax[i]);
251 }
252 #endif
253
254 BFT_MALLOC(xam, n_i_faces, cs_real_t);
255 BFT_MALLOC(dam, n_cells_ext, cs_real_33_t);
256 BFT_MALLOC(rhs, n_cells_ext, cs_real_3_t);
257
258 for (cs_lnum_t cell_id = 0; cell_id < n_cells_ext; cell_id++) {
259 for (int i = 0; i < 3; i++) {
260 for (int j = 0; j < 3; j++) {
261 dam[cell_id][i][j] = 0.;
262 }
263 rhs[cell_id][i] = 0.;
264 }
265 }
266
267 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
268 cs_lnum_t cell_id1 = i_face_cells[face_id][0];
269 cs_lnum_t cell_id2 = i_face_cells[face_id][1];
270 xam[face_id] = 0.;
271
272 //FIXME usefull?
273 double surf = surfn[face_id];
274 double vol = 0.5 * (volume[cell_id1] + volume[cell_id2]);
275 surf = CS_MAX(surf, 0.1*vol/dist[face_id]);
276 double ssd = surf / dist[face_id];
277
278 for (int i = 0; i < 3; i++) {
279 dam[cell_id1][i][i] += ssd;
280 dam[cell_id2][i][i] += ssd;
281 }
282
283 if ( mq->bad_cell_flag[cell_id1] & CS_BAD_CELL_TO_REGULARIZE
284 && mq->bad_cell_flag[cell_id2] & CS_BAD_CELL_TO_REGULARIZE) {
285 xam[face_id] = -ssd;
286 }
287 else if (mq->bad_cell_flag[cell_id1] & CS_BAD_CELL_TO_REGULARIZE) {
288 for (int i = 0; i < 3; i++) {
289 rhs[cell_id1][i] += ssd * var[cell_id2][i];
290 rhs[cell_id2][i] += ssd * var[cell_id2][i];
291 }
292 }
293 else if (mq->bad_cell_flag[cell_id2] & CS_BAD_CELL_TO_REGULARIZE) {
294 for (int i = 0; i < 3; i++) {
295 rhs[cell_id2][i] += ssd * var[cell_id1][i];
296 rhs[cell_id1][i] += ssd * var[cell_id1][i];
297 }
298 }
299 else {
300 for (int i = 0; i < 3; i++) {
301 rhs[cell_id1][i] += ssd * var[cell_id1][i];
302 rhs[cell_id2][i] += ssd * var[cell_id2][i];
303 }
304 }
305 }
306
307 /* Boudanry projection... should be consistent with BCs... */
308 if (boundary_projection == 1) {
309 for (cs_lnum_t face_id = 0; face_id < n_b_faces; face_id++) {
310 if (cs_glob_bc_type[face_id] == CS_SMOOTHWALL ||
311 cs_glob_bc_type[face_id] == CS_ROUGHWALL ||
312 cs_glob_bc_type[face_id] == CS_SYMMETRY ) {
313 cs_lnum_t cell_id = b_face_cells[face_id];
314 if (mq->bad_cell_flag[cell_id] & CS_BAD_CELL_TO_REGULARIZE) {
315 double ssd = surfbn[face_id] / distbr[face_id];
316 for (int i = 0; i < 3; i++) {
317 for (int j = 0; j < 3; j++) {
318 double nn = surfbo[face_id][i]/surfbn[face_id]
319 * surfbo[face_id][j]/surfbn[face_id];
320 dam[cell_id][i][j] += ssd * nn;
321 }
322 }
323 }
324 }
325 }
326 }
327
328 cs_real_t rnorm = sqrt(cs_gdot(3*n_cells,
329 (const cs_real_t *)rhs,
330 (const cs_real_t *)rhs));
331
332 /* Solver residual */
333 double ressol = 0.;
334 int niterf = 0;
335
336 cs_real_t epsilp = 1.e-12;
337
338 /* Matrix block size */
339 cs_lnum_t db_size[4] = {3, 3, 3, 3*3};
340
341 cs_sles_solve_native(-1, /* f_id */
342 "potential_regularisation_vector",
343 true, /* symmetric */
344 db_size,
345 NULL, /* eb_size */
346 (cs_real_t *)dam,
347 xam,
348 epsilp,
349 rnorm,
350 &niterf,
351 &ressol,
352 (cs_real_t *)rhs,
353 (cs_real_t *)var);
354
355 bft_printf("Solving %s: N iter: %d, Res: %12.5e, Norm: %12.5e\n",
356 "potential_regularisation_vector", niterf, ressol, rnorm);
357
358 //FIXME useless clipping? the matrix is min/max preserving..
359 #if 1
360 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
361 for (int i = 0; i < 3; i++) {
362 var[cell_id][i] = CS_MIN(var[cell_id][i], varmax[i]);
363 var[cell_id][i] = CS_MAX(var[cell_id][i], varmin[i]);
364 }
365 }
366 #endif
367
368 //FIXME periodicity of rotation
369 if (mesh->halo != NULL)
370 cs_halo_sync_var_strided(mesh->halo, CS_HALO_STANDARD, (cs_real_t *)var, 3);
371
372 /* Free solver setup */
373 cs_sles_free_native(-1, /* f_id*/
374 "potential_regularisation_vector");
375
376 /* Free memory */
377 BFT_FREE(xam);
378 BFT_FREE(dam);
379 BFT_FREE(rhs);
380 }
381
382 /*----------------------------------------------------------------------------*/
383 /*!
384 * \brief Regularisation on bad cells for symmetric tensors.
385 */
386 /*----------------------------------------------------------------------------*/
387
388 void
cs_bad_cells_regularisation_sym_tensor(cs_real_6_t * var,int boundary_projection)389 cs_bad_cells_regularisation_sym_tensor(cs_real_6_t *var,
390 int boundary_projection)
391 {
392 const cs_mesh_t *mesh = cs_glob_mesh;
393 cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
394
395 if (!(cs_glob_mesh_quantities_flag & CS_BAD_CELLS_REGULARISATION))
396 return;
397
398 cs_lnum_t n_cells_ext = mesh->n_cells_with_ghosts;
399 cs_lnum_t n_cells = mesh->n_cells;
400 cs_lnum_t n_i_faces = mesh->n_i_faces;
401 const cs_lnum_2_t *i_face_cells = (const cs_lnum_2_t *)mesh->i_face_cells;
402
403 const cs_real_t *surfn = mq->i_face_surf;
404 double *dist = mq->i_dist;
405 double *volume = mq->cell_vol;
406
407 cs_real_66_t *dam;
408 cs_real_6_t *rhs;
409 cs_real_t *xam;
410 #if 1
411 double varmin[6] = { 1.e20, 1.e20, 1.e20, 1.e20, 1.e20, 1.e20};
412 double varmax[6] = {-1.e20, -1.e20,-1.e20, -1.e20, -1.e20,-1.e20};
413
414 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++)
415 if (!(mq->bad_cell_flag[cell_id] & CS_BAD_CELL_TO_REGULARIZE)) {
416 for (int i = 0; i < 6; i++) {
417 varmin[i] = CS_MIN(varmin[i], var[cell_id][i]);
418 varmax[i] = CS_MAX(varmax[i], var[cell_id][i]);
419 }
420 }
421
422 for (int i = 0; i < 6; i++) {
423 cs_parall_min(1, CS_DOUBLE, &varmin[i]);
424 cs_parall_max(1, CS_DOUBLE, &varmax[i]);
425 }
426 #endif
427
428 BFT_MALLOC(xam, n_i_faces, cs_real_t);
429 BFT_MALLOC(dam, n_cells_ext, cs_real_66_t);
430 BFT_MALLOC(rhs, n_cells_ext, cs_real_6_t);
431
432 for (cs_lnum_t cell_id = 0; cell_id < n_cells_ext; cell_id++) {
433 for (int i = 0; i < 6; i++) {
434 for (int j = 0; j < 6; j++) {
435 dam[cell_id][i][j] = 0.;
436 }
437 rhs[cell_id][i] = 0.;
438 }
439 }
440
441 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
442 cs_lnum_t cell_id1 = i_face_cells[face_id][0];
443 cs_lnum_t cell_id2 = i_face_cells[face_id][1];
444 xam[face_id] = 0.;
445
446 //FIXME usefull?
447 double surf = surfn[face_id];
448 double vol = 0.5 * (volume[cell_id1] + volume[cell_id2]);
449 surf = CS_MAX(surf, 0.1*vol/dist[face_id]);
450 double ssd = surf / dist[face_id];
451
452 for (int i = 0; i < 6; i++) {
453 dam[cell_id1][i][i] += ssd;
454 dam[cell_id2][i][i] += ssd;
455 }
456
457 if ( mq->bad_cell_flag[cell_id1] & CS_BAD_CELL_TO_REGULARIZE
458 && mq->bad_cell_flag[cell_id2] & CS_BAD_CELL_TO_REGULARIZE) {
459 xam[face_id] = -ssd;
460 }
461 else if (mq->bad_cell_flag[cell_id1] & CS_BAD_CELL_TO_REGULARIZE) {
462 for (int i = 0; i < 6; i++) {
463 rhs[cell_id1][i] += ssd * var[cell_id2][i];
464 rhs[cell_id2][i] += ssd * var[cell_id2][i];
465 }
466 }
467 else if (mq->bad_cell_flag[cell_id2] & CS_BAD_CELL_TO_REGULARIZE) {
468 for (int i = 0; i < 6; i++) {
469 rhs[cell_id2][i] += ssd * var[cell_id1][i];
470 rhs[cell_id1][i] += ssd * var[cell_id1][i];
471 }
472 }
473 else {
474 for (int i = 0; i < 6; i++) {
475 rhs[cell_id1][i] += ssd * var[cell_id1][i];
476 rhs[cell_id2][i] += ssd * var[cell_id2][i];
477 }
478 }
479 }
480
481 /* Boundary projection... should be consistent with BCs... */
482 #if 0
483 if (boundary_projection == 1) {
484 cs_lnum_t n_b_faces = mesh->n_b_faces;
485 const cs_lnum_t *b_face_cells = mesh->b_face_cells;
486 const cs_real_3_t *surfbo = (const cs_real_3_t *) mq->b_face_normal;
487 const cs_real_t *surfbn = mq->b_face_surf;
488 double *distbr = mq->b_dist;
489
490 for (cs_lnum_t face_id = 0; face_id < n_b_faces; face_id++) {
491 if (cs_glob_bc_type[face_id] == CS_SMOOTHWALL ||
492 cs_glob_bc_type[face_id] == CS_ROUGHWALL ||
493 cs_glob_bc_type[face_id] == CS_SYMMETRY ) {
494 cs_lnum_t cell_id = b_face_cells[face_id];
495 if (mq->bad_cell_flag[cell_id] & CS_BAD_CELL_TO_REGULARIZE) {
496 double ssd = surfbn[face_id] / distbr[face_id];
497 for (int i = 0; i < 3; i++) {
498 for (int j = 0; j < 3; j++) {
499 double nn = surfbo[face_id][i]/surfbn[face_id]
500 * surfbo[face_id][j]/surfbn[face_id];
501 //TODO ??? dam[cell_id][i][j] += ssd * nn;
502 }
503 }
504 }
505 }
506 }
507 }
508 #endif
509
510 cs_real_t rnorm = sqrt(cs_gdot(6*n_cells,
511 (const cs_real_t *)rhs,
512 (const cs_real_t *)rhs));
513
514 /* Solver residual */
515 double ressol = 0.;
516 int niterf = 0;
517
518 cs_real_t epsilp = 1.e-12;
519
520 /* Matrix block size */
521 cs_lnum_t db_size[4] = {6, 6, 6, 6*6};
522
523 cs_sles_solve_native(-1, /* f_id */
524 "potential_regularisation_sym_tensor",
525 true, /* symmetric */
526 db_size,
527 NULL, /* eb_size */
528 (cs_real_t *)dam,
529 xam,
530 epsilp,
531 rnorm,
532 &niterf,
533 &ressol,
534 (cs_real_t *)rhs,
535 (cs_real_t *)var);
536
537 bft_printf("Solving %s: N iter: %d, Res: %12.5e, Norm: %12.5e\n",
538 "potential_regularisation_sym_tensor", niterf, ressol, rnorm);
539 //FIXME useless clipping? the matrix is min/max preserving..
540 #if 1
541 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
542 for (int i = 0; i < 6; i++) {
543 var[cell_id][i] = CS_MIN(var[cell_id][i], varmax[i]);
544 var[cell_id][i] = CS_MAX(var[cell_id][i], varmin[i]);
545 }
546 }
547 #endif
548
549 //FIXME periodicity of rotation
550 if (mesh->halo != NULL)
551 cs_halo_sync_var_strided(mesh->halo, CS_HALO_STANDARD, (cs_real_t *)var, 6);
552
553 /* Free solver setup */
554 cs_sles_free_native(-1, /* f_id*/
555 "potential_regularisation_sym_tensor");
556
557 /* Free memory */
558 BFT_FREE(xam);
559 BFT_FREE(dam);
560 BFT_FREE(rhs);
561 }
562
563 /*----------------------------------------------------------------------------*/
564 /*!
565 * \brief Regularisation on bad cells for tensors
566 *
567 */
568 /*----------------------------------------------------------------------------*/
569
570 void
cs_bad_cells_regularisation_tensor(cs_real_9_t * var,int boundary_projection)571 cs_bad_cells_regularisation_tensor(cs_real_9_t *var,
572 int boundary_projection)
573 {
574 const cs_mesh_t *mesh = cs_glob_mesh;
575 cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
576
577 if (!(cs_glob_mesh_quantities_flag & CS_BAD_CELLS_REGULARISATION))
578 return;
579
580 cs_lnum_t n_cells_ext = mesh->n_cells_with_ghosts;
581 cs_lnum_t n_cells = mesh->n_cells;
582 cs_lnum_t n_i_faces = mesh->n_i_faces;
583 const cs_lnum_2_t *i_face_cells = (const cs_lnum_2_t *)mesh->i_face_cells;
584
585 const cs_real_t *surfn = mq->i_face_surf;
586 double *dist = mq->i_dist;
587 double *volume = mq->cell_vol;
588
589 cs_real_99_t *dam;
590 cs_real_9_t *rhs;
591 cs_real_t *xam;
592 #if 1
593 double varmin[9] = { 1.e20, 1.e20, 1.e20, 1.e20, 1.e20, 1.e20, 1.e20, 1.e20, 1.e20};
594 double varmax[9] = {-1.e20, -1.e20,-1.e20, -1.e20, -1.e20,-1.e20, -1.e20, -1.e20,-1.e20};
595
596 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++)
597 if (!(mq->bad_cell_flag[cell_id] & CS_BAD_CELL_TO_REGULARIZE)) {
598 for (int i = 0; i < 9; i++) {
599 varmin[i] = CS_MIN(varmin[i], var[cell_id][i]);
600 varmax[i] = CS_MAX(varmax[i], var[cell_id][i]);
601 }
602 }
603
604 for (int i = 0; i < 9; i++) {
605 cs_parall_min(1, CS_DOUBLE, &varmin[i]);
606 cs_parall_max(1, CS_DOUBLE, &varmax[i]);
607 }
608 #endif
609
610 BFT_MALLOC(xam, n_i_faces, cs_real_t);
611 BFT_MALLOC(dam, n_cells_ext, cs_real_99_t);
612 BFT_MALLOC(rhs, n_cells_ext, cs_real_9_t);
613
614 for (cs_lnum_t cell_id = 0; cell_id < n_cells_ext; cell_id++) {
615 for (int i = 0; i < 9; i++) {
616 for (int j = 0; j < 9; j++) {
617 dam[cell_id][i][j] = 0.;
618 }
619 rhs[cell_id][i] = 0.;
620 }
621 }
622
623 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
624 cs_lnum_t cell_id1 = i_face_cells[face_id][0];
625 cs_lnum_t cell_id2 = i_face_cells[face_id][1];
626 xam[face_id] = 0.;
627
628 //FIXME usefull?
629 double surf = surfn[face_id];
630 double vol = 0.5 * (volume[cell_id1] + volume[cell_id2]);
631 surf = CS_MAX(surf, 0.1*vol/dist[face_id]);
632 double ssd = surf / dist[face_id];
633
634 for (int i = 0; i < 9; i++) {
635 dam[cell_id1][i][i] += ssd;
636 dam[cell_id2][i][i] += ssd;
637 }
638
639 if ( mq->bad_cell_flag[cell_id1] & CS_BAD_CELL_TO_REGULARIZE
640 && mq->bad_cell_flag[cell_id2] & CS_BAD_CELL_TO_REGULARIZE) {
641 xam[face_id] = -ssd;
642 }
643 else if (mq->bad_cell_flag[cell_id1] & CS_BAD_CELL_TO_REGULARIZE) {
644 for (int i = 0; i < 9; i++) {
645 rhs[cell_id1][i] += ssd * var[cell_id2][i];
646 rhs[cell_id2][i] += ssd * var[cell_id2][i];
647 }
648 }
649 else if (mq->bad_cell_flag[cell_id2] & CS_BAD_CELL_TO_REGULARIZE) {
650 for (int i = 0; i < 9; i++) {
651 rhs[cell_id2][i] += ssd * var[cell_id1][i];
652 rhs[cell_id1][i] += ssd * var[cell_id1][i];
653 }
654 }
655 else {
656 for (int i = 0; i < 9; i++) {
657 rhs[cell_id1][i] += ssd * var[cell_id1][i];
658 rhs[cell_id2][i] += ssd * var[cell_id2][i];
659 }
660 }
661 }
662
663 #if 0
664 /* Boudanry projection... should be consistent with BCs... */
665 if (boundary_projection == 1) {
666 const cs_lnum_t n_b_faces = mesh->n_b_faces;
667 const cs_real_t *distbr = mq->b_dist;
668 const cs_real_t *surfbn = mq->b_face_surf;
669 const cs_real_3_t *surfbo = (const cs_real_3_t *) mq->b_face_normal;
670 const cs_lnum_t *b_face_cells = mesh->b_face_cells;
671 for (cs_lnum_t face_id = 0; face_id < n_b_faces; face_id++) {
672 if (cs_glob_bc_type[face_id] == CS_SMOOTHWALL ||
673 cs_glob_bc_type[face_id] == CS_ROUGHWALL ||
674 cs_glob_bc_type[face_id] == CS_SYMMETRY ) {
675 cs_lnum_t cell_id = b_face_cells[face_id];
676 if (mq->bad_cell_flag[cell_id] & CS_BAD_CELL_TO_REGULARIZE) {
677 double ssd = surfbn[face_id] / distbr[face_id];
678 for (int i = 0; i < 3; i++) {
679 for (int j = 0; j < 3; j++) {
680 double nn = surfbo[face_id][i]/surfbn[face_id]
681 * surfbo[face_id][j]/surfbn[face_id];
682 //TODO ??? dam[cell_id][i][j] += ssd * nn;
683 }
684 }
685 }
686 }
687 }
688 }
689 #endif
690
691 cs_real_t rnorm = sqrt(cs_gdot(9*n_cells,
692 (const cs_real_t *)rhs,
693 (const cs_real_t *)rhs));
694
695 /* Solver residual */
696 double ressol = 0.;
697 int niterf = 0;
698
699 cs_real_t epsilp = 1.e-12;
700
701 /* Matrix block size */
702 cs_lnum_t db_size[4] = {9, 9, 9, 9*9};
703
704 cs_sles_solve_native(-1, /* f_id */
705 "potential_regularisation_tensor",
706 true, /* symmetric */
707 db_size,
708 NULL, /* eb_size */
709 (cs_real_t *)dam,
710 xam,
711 epsilp,
712 rnorm,
713 &niterf,
714 &ressol,
715 (cs_real_t *)rhs,
716 (cs_real_t *)var);
717
718 bft_printf("Solving %s: N iter: %d, Res: %12.5e, Norm: %12.5e\n",
719 "potential_regularisation_tensor", niterf, ressol, rnorm);
720 //FIXME useless clipping? the matrix is min/max preserving..
721 #if 1
722 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
723 for (int i = 0; i < 9; i++) {
724 var[cell_id][i] = CS_MIN(var[cell_id][i], varmax[i]);
725 var[cell_id][i] = CS_MAX(var[cell_id][i], varmin[i]);
726 }
727 }
728 #endif
729
730 //FIXME periodicity of rotation
731 if (mesh->halo != NULL)
732 cs_halo_sync_var_strided(mesh->halo, CS_HALO_STANDARD, (cs_real_t *)var, 9);
733
734 /* Free solver setup */
735 cs_sles_free_native(-1, /* f_id*/
736 "potential_regularisation_tensor");
737
738 /* Free memory */
739 BFT_FREE(xam);
740 BFT_FREE(dam);
741 BFT_FREE(rhs);
742 }
743
744 /*----------------------------------------------------------------------------*/
745
746 END_C_DECLS
747