1 /*============================================================================
2 * Log field and other array statistics at relevant time steps.
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 <math.h>
35 #include <stdarg.h>
36 #include <stdio.h>
37 #include <stdlib.h>
38 #include <string.h>
39
40 /*----------------------------------------------------------------------------
41 * Local headers
42 *----------------------------------------------------------------------------*/
43
44 #include "bft_mem.h"
45 #include "bft_printf.h"
46
47 #include "cs_array_reduce.h"
48 #include "cs_base.h"
49 #include "cs_ctwr.h"
50 #include "cs_fan.h"
51 #include "cs_field.h"
52 #include "cs_log.h"
53 #include "cs_map.h"
54 #include "cs_mesh.h"
55 #include "cs_mesh_location.h"
56 #include "cs_parall.h"
57 #include "cs_parameters.h"
58 #include "cs_prototypes.h"
59 #include "cs_range_set.h"
60 #include "cs_time_moment.h"
61 #include "cs_time_plot.h"
62 #include "cs_time_step.h"
63 #include "cs_lagr_stat.h"
64 #include "cs_lagr_log.h"
65
66 /*----------------------------------------------------------------------------
67 * Header for the current file
68 *----------------------------------------------------------------------------*/
69
70 #include "cs_log_iteration.h"
71
72 /*----------------------------------------------------------------------------*/
73
74 BEGIN_C_DECLS
75
76 /*=============================================================================
77 * Additional doxygen documentation
78 *============================================================================*/
79
80 /*!
81 \file cs_log_iteration.c
82
83 \brief Log field and other array statistics at relevant time steps.
84 */
85
86 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
87
88 /*============================================================================
89 * Type and macro definitions
90 *============================================================================*/
91
92 /* Simple statistics */
93 /*-------------------*/
94
95 typedef struct {
96
97 int name_id; /* Associated name id */
98 int cat_id; /* Associated category id */
99 int loc_id; /* Associated mesh location id */
100 bool intensive; /* Are associated values intensive ? */
101 int dim; /* Associated dimension */
102 int v_idx; /* Start index of values */
103
104 } cs_log_sstats_t;
105
106 /* Clipping info */
107 /*---------------*/
108
109 typedef struct {
110
111 int f_id; /* Associated field id, or -1 */
112 int name_id; /* Associated name id if not a field, -1 for a field */
113 int dim; /* Associated dimension */
114 int c_idx; /* Start index of counts */
115 int v_idx; /* Start index of values */
116
117 } cs_log_clip_t;
118
119 /*============================================================================
120 * Static global variables
121 *============================================================================*/
122
123 static cs_map_name_to_id_t *_name_map = NULL;
124
125 static cs_map_name_to_id_t *_category_map = NULL;
126 static int _sstats_val_size = 0;
127 static int _sstats_val_size_max = 0;
128 static int _n_sstats = 0;
129 static int _n_sstats_max = 0;
130 static double *_sstats_vmin = NULL;
131 static double *_sstats_vmax = NULL;
132 static double *_sstats_vsum = NULL;
133 static double *_sstats_wsum = NULL;
134 static cs_log_sstats_t *_sstats = NULL;
135
136 static int _clips_val_size = 0;
137 static int _clips_val_size_max = 0;
138 static int _n_clips = 0;
139 static int _n_clips_max = 0;
140 static cs_gnum_t *_clips_count = NULL;
141 static double *_clips_vmin = NULL;
142 static double *_clips_vmax = NULL;
143 static cs_log_clip_t *_clips = NULL;
144
145 static cs_time_plot_t *_l2_residual_plot = NULL;
146
147 /*============================================================================
148 * Prototypes for functions intended for use only by Fortran wrappers.
149 * (descriptions follow, with function bodies).
150 *============================================================================*/
151
152 /*============================================================================
153 * Fortran function prototypes for subroutines from field.f90.
154 *============================================================================*/
155
156 /*============================================================================
157 * Private function definitions
158 *============================================================================*/
159
160 /*----------------------------------------------------------------------------
161 * Compare simple stats elements (qsort function).
162 *
163 * parameters:
164 * x <-> pointer to first element
165 * y <-> pointer to second element
166 *
167 * returns:
168 * -1 if x < y, 0 if x = y, or 1 if x > y
169 *----------------------------------------------------------------------------*/
170
171 static int
_compare_sstats(const void * x,const void * y)172 _compare_sstats(const void *x,
173 const void *y)
174 {
175 int retval = 1;
176
177 const cs_log_sstats_t *s0 = x;
178 const cs_log_sstats_t *s1 = y;
179
180 if (s0->cat_id < s1->cat_id)
181 retval = -1;
182
183 else if (s0->cat_id == s1->cat_id) {
184 if (s0->name_id < s1->name_id)
185 retval = -1;
186 else if (s0->name_id == s1->name_id)
187 retval = 0;
188 }
189
190 return retval;
191 }
192
193 /*----------------------------------------------------------------------------
194 * Find simple stats by id
195 *
196 * parameters:
197 * cat_id <-- category id
198 * name_id <-- name id
199 *
200 * returns:
201 * id of simple stats, or -1 if not found
202 *----------------------------------------------------------------------------*/
203
204 static int
_find_sstats(int cat_id,int name_id)205 _find_sstats(int cat_id,
206 int name_id)
207 {
208 /* Use binary search to find entry */
209
210 int sstat_id = -1;
211
212 if (_n_sstats > 0) {
213
214 int start_id = 0;
215 int end_id = _n_sstats;
216
217 while (start_id < end_id) {
218 int mid_id = start_id + ((end_id -start_id) / 2);
219 int cmp_ret = 0;
220 int cmp_cat = _sstats[mid_id].cat_id;
221 if (cmp_cat < cat_id)
222 cmp_ret = -1;
223 else if (cmp_cat > cat_id)
224 cmp_ret = 1;
225 else {
226 int cmp_name = _sstats[mid_id].name_id;
227 if (cmp_name < name_id)
228 cmp_ret = -1;
229 else if (cmp_name > name_id)
230 cmp_ret = 1;
231 else {
232 sstat_id = mid_id;
233 break;
234 }
235 }
236 if (cmp_ret < 0)
237 start_id = mid_id + 1;
238 else if (cmp_ret > 0)
239 end_id = mid_id;
240 }
241 }
242
243 return sstat_id;
244 }
245
246 /*----------------------------------------------------------------------------
247 * Compare clipping elements (qsort function).
248 *
249 * parameters:
250 * x <-> pointer to first element
251 * y <-> pointer to second element
252 *
253 * returns:
254 * -1 if x < y, 0 if x = y, or 1 if x > y
255 *----------------------------------------------------------------------------*/
256
_compare_clips(const void * x,const void * y)257 static int _compare_clips(const void *x, const void *y)
258 {
259 int retval = 1;
260
261 const cs_log_clip_t *c0 = x;
262 const cs_log_clip_t *c1 = y;
263
264 if (c0->name_id < c1->name_id)
265 retval = -1;
266
267 else if (c0->name_id == c1->name_id) {
268 if (c0->f_id < c1->f_id)
269 retval = -1;
270 else if (c0->f_id == c1->f_id)
271 retval = 0;
272 }
273
274 return retval;
275 }
276
277 /*----------------------------------------------------------------------------
278 * Find clippings by id
279 *
280 * parameters:
281 * f_id <-- field id
282 * name_id <-- name id
283 *
284 * returns:
285 * id of clipping, or -1 if not found
286 *----------------------------------------------------------------------------*/
287
288 static int
_find_clip(int f_id,int name_id)289 _find_clip(int f_id,
290 int name_id)
291 {
292 /* Use binary search to find entry */
293
294 int clip_id = -1;
295
296 if (_n_clips > 0) {
297
298 int start_id = 0;
299 int end_id = _n_clips;
300
301 while (start_id < end_id) {
302 int mid_id = start_id + ((end_id -start_id) / 2);
303 int cmp_ret = 0;
304 int cmp_name = _clips[mid_id].name_id;
305 if (cmp_name < name_id)
306 cmp_ret = -1;
307 else if (cmp_name > name_id)
308 cmp_ret = 1;
309 else {
310 int cmp_f = _clips[mid_id].f_id;
311 if (cmp_f < f_id)
312 cmp_ret = -1;
313 else if (cmp_f > f_id)
314 cmp_ret = 1;
315 else {
316 clip_id = mid_id;
317 break;
318 }
319 }
320 if (cmp_ret < 0)
321 start_id = mid_id + 1;
322 else if (cmp_ret > 0)
323 end_id = mid_id;
324 }
325 }
326
327 return clip_id;
328 }
329
330 /*----------------------------------------------------------------------------
331 * Log information for a given array.
332 *
333 * parameters:
334 * prefix <-- string inserted before name
335 * name <-- array name or label
336 * name_width <-- width of "name" column
337 * dim <-- array dimension
338 * n_g_elts <-- global number of associated elements,
339 * total_weight <-- if > 0, weight (volume or surface) of array location
340 * vmin <-- minimum values of each component or norm
341 * vmax <-- maximum values of each component or norm
342 * vsum <-- sum of each component or norm
343 * wsum <-- weighted sum of each component or norm, or NULL
344 * fpe_flag <-- was a "not a number" or floating-point error detected ?
345 *----------------------------------------------------------------------------*/
346
347 static void
_log_array_info(const char * prefix,const char * name,size_t name_width,int dim,int n_g_elts,double total_weight,double vmin[],const double vmax[],const double vsum[],const double * wsum,int * fpe_flag)348 _log_array_info(const char *prefix,
349 const char *name,
350 size_t name_width,
351 int dim,
352 int n_g_elts,
353 double total_weight,
354 double vmin[],
355 const double vmax[],
356 const double vsum[],
357 const double *wsum,
358 int *fpe_flag)
359 {
360 const int _dim = (dim == 3) ? 4 : dim;
361
362 char tmp_s[2][64] = {"", ""};
363
364 for (int c_id = 0; c_id < _dim; c_id++) {
365
366 if (dim == 3) {
367 if (c_id < 3)
368 snprintf(tmp_s[1], 63, "%s%s", name, cs_glob_field_comp_name_3[c_id]);
369 else
370 snprintf(tmp_s[1], 63, "ǁ%sǁ", name);
371 tmp_s[1][63] = '\0';
372 cs_log_strpad(tmp_s[0], tmp_s[1], name_width, 64);
373 }
374 else if (dim == 6) {
375 snprintf(tmp_s[1], 63, "%s%s", name, cs_glob_field_comp_name_6[c_id]);
376 tmp_s[1][63] = '\0';
377 cs_log_strpad(tmp_s[0], tmp_s[1], name_width, 64);
378 }
379 else if (dim == 9) {
380 snprintf(tmp_s[1], 63, "%s%s", name, cs_glob_field_comp_name_9[c_id]);
381 tmp_s[1][63] = '\0';
382 cs_log_strpad(tmp_s[0], tmp_s[1], name_width, 64);
383 }
384 else
385 cs_log_strpad(tmp_s[0], name, name_width, 64);
386
387 if (total_weight >= 0)
388 cs_log_printf(CS_LOG_DEFAULT,
389 "%s%s %14.5g %14.5g %14.5g %14.5g\n",
390 prefix,
391 tmp_s[0],
392 vmin[c_id],
393 vmax[c_id],
394 vsum[c_id] / n_g_elts,
395 wsum[c_id] / total_weight);
396 else
397 cs_log_printf(CS_LOG_DEFAULT,
398 "%s%s %14.5g %14.5g %14.5g\n",
399 prefix,
400 tmp_s[0],
401 vmin[c_id],
402 vmax[c_id],
403 vsum[c_id] / n_g_elts);
404
405 /* Check NAN and exit */
406 if (isnan(vsum[c_id]))
407 *fpe_flag = 1;
408 }
409
410 }
411
412 /*----------------------------------------------------------------------------
413 * Log information for a given clipping.
414 *
415 * parameters:
416 * prefix <-- string inserted before name
417 * name <-- array name or label
418 * name_width <-- width of "name" column
419 * dim <-- array dimension
420 * count_min <-- global number of clips to minimum
421 * count_max <-- global number of clips to maximum
422 * vmin <-- minimum values of each component
423 * vmax <-- maximum values of each component
424 *----------------------------------------------------------------------------*/
425
426 static void
_log_clip_info(const char * prefix,const char * name,size_t name_width,int dim,cs_gnum_t count_min[],cs_gnum_t count_max[],const double vmin[],const double vmax[])427 _log_clip_info(const char *prefix,
428 const char *name,
429 size_t name_width,
430 int dim,
431 cs_gnum_t count_min[],
432 cs_gnum_t count_max[],
433 const double vmin[],
434 const double vmax[])
435 {
436 int c_id;
437 const int _dim = (dim == 1) ? dim : dim + 1;
438
439 char tmp_s[2][64] = {"", ""};
440
441 for (c_id = 0; c_id < _dim; c_id++) {
442
443 int _count_max, _count_min;
444 double _vmin, _vmax;
445 const char *_name = (c_id == 0) ? name : " ";
446
447 if (c_id < _dim) {
448 _count_min = count_min[2*c_id];
449 _count_max = count_max[2*c_id];
450 _vmin = vmin[c_id];
451 _vmax = vmax[c_id];
452 }
453
454 if (dim > 1) {
455 if (c_id == 0) {
456 snprintf(tmp_s[1], 63, "%s", _name);
457 tmp_s[1][63] = '\0';
458 cs_log_strpad(tmp_s[0], tmp_s[1], name_width, 64);
459 _vmin = 0.;
460 _vmax = 0.;
461 for (int i = 0; i< dim; i++) {
462 _vmin += vmin[i]*vmin[i];
463 _vmax += vmax[i]*vmax[i];
464 }
465 _vmin = sqrt(_vmin);
466 _vmax = sqrt(_vmax);
467 }
468 else {
469 if (dim == 3)
470 snprintf(tmp_s[1], 63, "%s%s", name, cs_glob_field_comp_name_3[c_id - 1]);
471 else if (dim == 6)
472 snprintf(tmp_s[1], 63, "%s%s", name, cs_glob_field_comp_name_6[c_id - 1]);
473 tmp_s[1][63] = '\0';
474 cs_log_strpad(tmp_s[0], tmp_s[1], name_width, 64);
475 }
476 }
477 else
478 cs_log_strpad(tmp_s[0], name, name_width, 64);
479
480 if (_count_min > 0 && _count_max > 0)
481 cs_log_printf(CS_LOG_DEFAULT,
482 "%s%s %14.5g %14.5g %12llu %12llu\n",
483 prefix,
484 tmp_s[0],
485 _vmin,
486 _vmax,
487 (unsigned long long)_count_min,
488 (unsigned long long)_count_max);
489 else if (_count_min > 0)
490 cs_log_printf(CS_LOG_DEFAULT,
491 "%s%s %14.5g %12llu\n",
492 prefix,
493 tmp_s[0],
494 _vmin,
495 (unsigned long long)_count_min);
496 else if (_count_max > 0)
497 cs_log_printf(CS_LOG_DEFAULT,
498 "%s%s %14.5g %12llu\n",
499 prefix,
500 tmp_s[0],
501 _vmax,
502 (unsigned long long)_count_max);
503 else
504 cs_log_printf(CS_LOG_DEFAULT,
505 "%s%s\n",
506 prefix,
507 tmp_s[0]);
508 }
509 }
510
511 /*----------------------------------------------------------------------------
512 * Main logging output of variables.
513 *----------------------------------------------------------------------------*/
514
515 static void
_log_fields(void)516 _log_fields(void)
517 {
518 int f_id, li, log_count;
519
520 int log_count_max = 0;
521 int fpe_flag = 0;
522 int *log_id = NULL, *moment_id = NULL;
523 double *vmin = NULL, *vmax = NULL, *vsum = NULL, *wsum = NULL;
524
525 char tmp_s[5][64] = {"", "", "", "", ""};
526
527 const char _underline[] = "---------------------------------";
528 const int n_fields = cs_field_n_fields();
529 const int n_moments = cs_time_moment_n_moments();
530 const int log_key_id = cs_field_key_id("log");
531 const int label_key_id = cs_field_key_id("label");
532
533 const cs_mesh_location_type_t m_l[] = {CS_MESH_LOCATION_CELLS,
534 CS_MESH_LOCATION_INTERIOR_FACES,
535 CS_MESH_LOCATION_BOUNDARY_FACES,
536 CS_MESH_LOCATION_VERTICES
537 };
538
539 cs_mesh_t *m = cs_glob_mesh;
540 const cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
541
542 /* Allocate working arrays */
543
544 log_count_max = n_fields;
545
546 BFT_MALLOC(log_id, n_fields, int);
547 BFT_MALLOC(vmin, log_count_max, double);
548 BFT_MALLOC(vmax, log_count_max, double);
549 BFT_MALLOC(vsum, log_count_max, double);
550 BFT_MALLOC(wsum, log_count_max, double);
551
552 if (n_moments > 0) {
553 BFT_MALLOC(moment_id, n_fields, int);
554 for (f_id = 0; f_id < n_fields; f_id++)
555 moment_id[f_id] = -1;
556 for (int m_id = 0; m_id < n_moments; m_id++) {
557 const cs_field_t *f = cs_time_moment_get_field(m_id);
558 if (f != NULL)
559 moment_id[f->id] = m_id;
560 }
561 }
562
563 /* Loop on locations */
564
565 for (li = 0; li < 4; li++) {
566
567 size_t max_name_width = cs_log_strlen(_("field"));
568 int loc_id = m_l[li];
569 cs_lnum_t have_weight = 0;
570 double total_weight = -1;
571 cs_gnum_t n_g_elts = 0;
572 cs_real_t *gather_array = NULL; /* only if CS_MESH_LOCATION_VERTICES */
573 const cs_lnum_t *n_elts = cs_mesh_location_get_n_elts(loc_id);
574 const cs_lnum_t _n_elts = n_elts[0];
575 const cs_real_t *weight = NULL;
576
577 if (mq != NULL) {
578 switch(loc_id) {
579 case CS_MESH_LOCATION_CELLS:
580 n_g_elts = m->n_g_cells;
581 weight = mq->cell_vol;
582 have_weight = 1;
583 total_weight = mq->tot_vol;
584 break;
585 case CS_MESH_LOCATION_INTERIOR_FACES:
586 n_g_elts = m->n_g_i_faces;
587 weight = mq->i_face_surf;
588 cs_array_reduce_sum_l(_n_elts, 1, NULL, weight, &total_weight);
589 cs_parall_sum(1, CS_DOUBLE, &total_weight);
590 have_weight = 1;
591 break;
592 case CS_MESH_LOCATION_BOUNDARY_FACES:
593 n_g_elts = m->n_g_b_faces;
594 weight = mq->b_face_surf;
595 cs_array_reduce_sum_l(_n_elts, 1, NULL, weight, &total_weight);
596 cs_parall_sum(1, CS_DOUBLE, &total_weight);
597 have_weight = 1;
598 break;
599 case CS_MESH_LOCATION_VERTICES:
600 n_g_elts = m->n_g_vertices;
601 have_weight = 0;
602 BFT_MALLOC(gather_array, m->n_vertices, cs_real_t);
603 break;
604 default:
605 n_g_elts = _n_elts;
606 cs_parall_counter(&n_g_elts, 1);
607 break;
608 }
609 }
610
611 if (n_g_elts == 0)
612 continue;
613
614 /* First loop on fields */
615
616 log_count = 0;
617
618 for (f_id = 0; f_id < n_fields; f_id++) {
619
620 int _dim, c_id;
621
622 const cs_field_t *f = cs_field_by_id(f_id);
623
624 if (f->location_id != loc_id || ! (cs_field_get_key_int(f, log_key_id))) {
625 log_id[f_id] = -1;
626 continue;
627 }
628
629 /* Only log active moments */
630
631 if (moment_id != NULL) {
632 if (moment_id[f_id] > -1) {
633 if (!cs_time_moment_is_active(moment_id[f_id])) {
634 log_id[f_id] = -1;
635 continue;
636 }
637 }
638 }
639
640 /* Position in log */
641
642 log_id[f_id] = log_count;
643
644 _dim = (f->dim == 3) ? 4 : f->dim;
645
646 while (log_count + _dim > log_count_max) {
647 log_count_max *= 2;
648 BFT_REALLOC(vmin, log_count_max, double);
649 BFT_REALLOC(vmax, log_count_max, double);
650 BFT_REALLOC(vsum, log_count_max, double);
651 BFT_REALLOC(wsum, log_count_max, double);
652 }
653
654 if (have_weight && (f->type & CS_FIELD_INTENSIVE)) {
655 cs_array_reduce_simple_stats_l_w(_n_elts,
656 f->dim,
657 NULL,
658 NULL,
659 f->val,
660 weight,
661 vmin + log_count,
662 vmax + log_count,
663 vsum + log_count,
664 wsum + log_count);
665
666 }
667 else {
668
669 cs_real_t *field_val = f->val;
670 cs_lnum_t _n_cur_elts = _n_elts;
671 if (gather_array != NULL) { /* Eliminate shared values whose local
672 rank is not owner and compact */
673
674 if (m->vtx_range_set == NULL)
675 m->vtx_range_set = cs_range_set_create(m->vtx_interfaces,
676 NULL,
677 m->n_vertices,
678 false, /* balance */
679 2, /* tr_ignore */
680 0); /* g_id_base */
681
682 if (f->dim > 1)
683 BFT_REALLOC(gather_array, (f->dim * m->n_vertices), cs_real_t);
684
685 cs_range_set_gather(m->vtx_range_set,
686 CS_REAL_TYPE,
687 f->dim,
688 f->val,
689 gather_array);
690 field_val = gather_array;
691 _n_cur_elts = m->vtx_range_set->n_elts[0];
692 }
693 cs_array_reduce_simple_stats_l(_n_cur_elts,
694 f->dim,
695 NULL,
696 field_val,
697 vmin + log_count,
698 vmax + log_count,
699 vsum + log_count);
700
701 if (have_weight) {
702 for (c_id = 0; c_id < _dim; c_id++)
703 wsum[log_count + c_id] = 0.;
704 }
705 }
706
707 log_count += _dim;
708
709 const char *name = cs_field_get_key_str(f, label_key_id);
710 if (name == NULL)
711 name = f->name;
712
713 size_t l_name_width = cs_log_strlen(name);
714 if (f->dim == 3)
715 l_name_width += 3;
716 else if (f->dim > 3)
717 l_name_width += 4;
718
719 max_name_width = CS_MAX(max_name_width, l_name_width);
720
721 } /* End of first loop on fields */
722
723 if (gather_array != NULL)
724 BFT_FREE(gather_array);
725
726 if (log_count < 1)
727 continue;
728
729 /* Group MPI operations if required */
730
731 cs_parall_min(log_count, CS_DOUBLE, vmin);
732 cs_parall_max(log_count, CS_DOUBLE, vmax);
733 cs_parall_sum(log_count, CS_DOUBLE, vsum);
734
735 # if defined(HAVE_MPI)
736 cs_parall_counter_max(&have_weight, 1);
737 if (have_weight)
738 cs_parall_sum(log_count, CS_DOUBLE, wsum);
739 # endif
740
741 /* Print headers */
742
743 max_name_width = CS_MIN(max_name_width, 63);
744
745 const char *loc_name = _(cs_mesh_location_get_name(loc_id));
746 size_t loc_name_w = cs_log_strlen(loc_name);
747
748 cs_log_printf(CS_LOG_DEFAULT,
749 _("\n"
750 " ** Field values on %s\n"
751 " ----------------%.*s\n"),
752 loc_name, (int)loc_name_w, _underline);
753
754 cs_log_strpad(tmp_s[0], _("field"), max_name_width, 64);
755 cs_log_strpadl(tmp_s[1], _("minimum"), 14, 64);
756 cs_log_strpadl(tmp_s[2], _("maximum"), 14, 64);
757 cs_log_strpadl(tmp_s[3], _("set mean"), 14, 64);
758 if (have_weight) {
759 cs_log_strpadl(tmp_s[4], _("spatial mean"), 14, 64);
760 cs_log_printf(CS_LOG_DEFAULT,
761 "\n %s %s %s %s %s\n",
762 tmp_s[0], tmp_s[1], tmp_s[2], tmp_s[3], tmp_s[4]);
763 }
764 else
765 cs_log_printf(CS_LOG_DEFAULT,
766 "\n %s %s %s %s\n",
767 tmp_s[0], tmp_s[1], tmp_s[2], tmp_s[3]);
768
769 /* Underline */
770
771 size_t n_cols = (have_weight) ? 5 : 4;
772
773 for (size_t col = 0; col < n_cols; col++) {
774 size_t i;
775 size_t w0 = (col == 0) ? max_name_width : 14;
776 for (i = 0; i < w0; i++)
777 tmp_s[col][i] = '-';
778 tmp_s[col][w0] = '\0';
779 }
780 if (have_weight) {
781 cs_log_printf(CS_LOG_DEFAULT,
782 "- %s %s %s %s %s\n",
783 tmp_s[0], tmp_s[1], tmp_s[2], tmp_s[3], tmp_s[4]);
784 }
785 else
786 cs_log_printf(CS_LOG_DEFAULT,
787 "- %s %s %s %s\n",
788 tmp_s[0], tmp_s[1], tmp_s[2], tmp_s[3]);
789
790 /* Second loop on fields */
791
792 log_count = 0;
793
794 for (f_id = 0; f_id < n_fields; f_id++) {
795
796 int _dim;
797 if (log_id[f_id] < 0)
798 continue;
799
800 const cs_field_t *f = cs_field_by_id(f_id);
801
802 /* Position in log */
803
804 _dim = (f->dim == 3) ? 4 : f->dim;
805
806 const char *name = cs_field_get_key_str(f, label_key_id);
807 if (name == NULL)
808 name = f->name;
809
810 double t_weight = -1;
811 if (total_weight > 0 && (f->type & CS_FIELD_INTENSIVE))
812 t_weight = total_weight;
813
814 char prefix[] = "v ";
815 if (moment_id != NULL) {
816 if (moment_id[f_id] > -1)
817 prefix[0] = 'm';
818 }
819 if (f->type & CS_FIELD_ACCUMULATOR)
820 prefix[0] = 'm';
821
822 _log_array_info(prefix,
823 name,
824 max_name_width,
825 f->dim,
826 n_g_elts,
827 t_weight,
828 vmin + log_count,
829 vmax + log_count,
830 vsum + log_count,
831 wsum + log_count,
832 &fpe_flag);
833
834 log_count += _dim;
835
836 } /* End of loop on fields */
837
838 } /* End of loop on mesh locations */
839
840 /* Check NaN and exit */
841 if (fpe_flag == 1)
842 bft_error(__FILE__, __LINE__, 0,
843 _("Invalid (not-a-number) values detected for a field."));
844
845 BFT_FREE(moment_id);
846 BFT_FREE(wsum);
847 BFT_FREE(vsum);
848 BFT_FREE(vmax);
849 BFT_FREE(vmin);
850 BFT_FREE(log_id);
851
852 cs_log_printf(CS_LOG_DEFAULT, "\n");
853 }
854
855 /*----------------------------------------------------------------------------
856 * Main logging output of additional simple statistics
857 *----------------------------------------------------------------------------*/
858
859 static void
_log_sstats(void)860 _log_sstats(void)
861 {
862 int stat_id;
863 int fpe_flag = 0;
864 double _boundary_surf = -1;
865 double _interior_surf = -1;
866 double *vmin = NULL, *vmax = NULL, *vsum = NULL, *wsum = NULL;
867
868 char tmp_s[5][64] = {"", "", "", "", ""};
869
870 const char _underline[] = "---------------------------------";
871
872 const cs_mesh_t *m = cs_glob_mesh;
873 const cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
874
875 /* Allocate working arrays */
876
877 BFT_MALLOC(vmin, _sstats_val_size, double);
878 BFT_MALLOC(vmax, _sstats_val_size, double);
879 BFT_MALLOC(vsum, _sstats_val_size, double);
880 BFT_MALLOC(wsum, _sstats_val_size, double);
881
882 memcpy(vmin, _sstats_vmin, _sstats_val_size*sizeof(double));
883 memcpy(vmax, _sstats_vmax, _sstats_val_size*sizeof(double));
884 memcpy(vsum, _sstats_vsum, _sstats_val_size*sizeof(double));
885 memcpy(wsum, _sstats_wsum, _sstats_val_size*sizeof(double));
886
887 /* Group MPI operations if required */
888
889 cs_parall_min(_sstats_val_size, CS_DOUBLE, vmin);
890 cs_parall_max(_sstats_val_size, CS_DOUBLE, vmax);
891 cs_parall_sum(_sstats_val_size, CS_DOUBLE, vsum);
892 cs_parall_sum(_sstats_val_size, CS_DOUBLE, wsum);
893
894 /* Loop on statistics */
895
896 int sstat_cat_start = 0;
897
898 while (sstat_cat_start < _n_sstats) {
899
900 int sstat_cat_end = sstat_cat_start;
901
902 while (sstat_cat_end < _n_sstats) {
903 if (_sstats[sstat_cat_end].cat_id != _sstats[sstat_cat_start].cat_id)
904 break;
905 else
906 sstat_cat_end ++;
907 }
908
909 const char *cat_name
910 = cs_map_name_to_id_reverse(_category_map,
911 _sstats[sstat_cat_start].cat_id);
912
913 const size_t cat_name_w = cs_log_strlen(_(cat_name));
914 size_t max_name_width = cat_name_w;
915
916 /* Now separate by mesh location */
917
918 int loc_min = cs_mesh_location_n_locations() + 1;
919 int loc_max = -1;
920
921 for (stat_id = sstat_cat_start; stat_id < sstat_cat_end; stat_id++) {
922 int loc_id = _sstats[stat_id].loc_id;
923 loc_min = CS_MIN(loc_min, loc_id);
924 loc_max = CS_MAX(loc_max, loc_id);
925 }
926
927 for (int loc_id = loc_min; loc_id <= loc_max; loc_id++) {
928
929 int n_loc_stats = 0;
930 for (stat_id = sstat_cat_start; stat_id < sstat_cat_end; stat_id++) {
931 if (_sstats[stat_id].loc_id == loc_id)
932 n_loc_stats++;
933 }
934 if (n_loc_stats == 0)
935 continue;
936
937 cs_gnum_t n_g_elts = 0;
938 int have_weight = 0;
939 double total_weight = -1;
940 const cs_lnum_t *n_elts = cs_mesh_location_get_n_elts(loc_id);;
941 const cs_lnum_t _n_elts = n_elts[0];
942 const cs_real_t *weight = NULL;
943 const char *loc_name = _(cs_mesh_location_get_name(loc_id));
944 size_t loc_name_w = cs_log_strlen(loc_name);
945
946 if (mq != NULL) {
947 switch(loc_id) {
948 case CS_MESH_LOCATION_CELLS:
949 n_g_elts = m->n_g_cells;
950 weight = mq->cell_vol;
951 have_weight = 1;
952 total_weight = mq->tot_vol;
953 break;
954 case CS_MESH_LOCATION_INTERIOR_FACES:
955 n_g_elts = m->n_g_i_faces;
956 weight = mq->i_face_surf;
957 if (_interior_surf < 0) {
958 cs_array_reduce_sum_l(_n_elts, 1, NULL, weight, &_interior_surf);
959 cs_parall_sum(1, CS_DOUBLE, &_interior_surf);
960 if (_interior_surf < 0) _interior_surf = 0; /* just to be safe */
961 }
962 total_weight = _interior_surf;
963 have_weight = 1;
964 break;
965 case CS_MESH_LOCATION_BOUNDARY_FACES:
966 n_g_elts = m->n_g_b_faces;
967 weight = mq->b_face_surf;
968 if (_boundary_surf < 0) {
969 cs_array_reduce_sum_l(_n_elts, 1, NULL, weight, &_boundary_surf);
970 cs_parall_sum(1, CS_DOUBLE, &_boundary_surf);
971 if (_boundary_surf < 0) _boundary_surf = 0; /* just to be safe */
972 }
973 total_weight = _boundary_surf;
974 have_weight = 1;
975 break;
976 case CS_MESH_LOCATION_VERTICES:
977 n_g_elts = m->n_g_vertices;
978 have_weight = 0;
979 break;
980 default:
981 n_g_elts = _n_elts;
982 cs_parall_counter(&n_g_elts, 1);
983 break;
984 }
985 }
986
987 for (stat_id = sstat_cat_start; stat_id < sstat_cat_end; stat_id++) {
988 if (_sstats[stat_id].loc_id == loc_id) {
989 const char *stat_name
990 = cs_map_name_to_id_reverse(_name_map, _sstats[stat_id].name_id);
991 size_t name_width = strlen(stat_name);
992 max_name_width = CS_MAX(max_name_width, name_width);
993 }
994 }
995 max_name_width = CS_MIN(max_name_width, 63);
996
997 cs_log_printf(CS_LOG_DEFAULT,
998 _("\n"
999 " ** Computed values on %s\n"
1000 " -------------------%.*s\n"),
1001 loc_name, (int)loc_name_w, _underline);
1002
1003 cs_log_strpad(tmp_s[0], _(cat_name), max_name_width, 64);
1004 cs_log_strpadl(tmp_s[1], _("minimum"), 14, 64);
1005 cs_log_strpadl(tmp_s[2], _("maximum"), 14, 64);
1006 cs_log_strpadl(tmp_s[3], _("set mean"), 14, 64);
1007 if (have_weight) {
1008 cs_log_strpadl(tmp_s[4], _("spatial mean"), 14, 64);
1009 cs_log_printf(CS_LOG_DEFAULT,
1010 "\n %s %s %s %s %s\n",
1011 tmp_s[0], tmp_s[1], tmp_s[2], tmp_s[3], tmp_s[4]);
1012 }
1013 else
1014 cs_log_printf(CS_LOG_DEFAULT,
1015 "\n %s %s %s %s\n",
1016 tmp_s[0], tmp_s[1], tmp_s[2], tmp_s[3]);
1017
1018 /* Underline */
1019
1020 size_t n_cols = (have_weight) ? 5 : 4;
1021
1022 for (size_t col = 0; col < n_cols; col++) {
1023 size_t i;
1024 size_t w0 = (col == 0) ? max_name_width : 14;
1025 for (i = 0; i < w0; i++)
1026 tmp_s[col][i] = '-';
1027 tmp_s[col][w0] = '\0';
1028 }
1029 if (have_weight) {
1030 cs_log_printf(CS_LOG_DEFAULT,
1031 " %s %s %s %s %s\n",
1032 tmp_s[0], tmp_s[1], tmp_s[2], tmp_s[3], tmp_s[4]);
1033 }
1034 else
1035 cs_log_printf(CS_LOG_DEFAULT,
1036 " %s %s %s %s\n",
1037 tmp_s[0], tmp_s[1], tmp_s[2], tmp_s[3]);
1038
1039 /* Print values */
1040
1041 for (stat_id = sstat_cat_start; stat_id < sstat_cat_end; stat_id++) {
1042
1043 if (_sstats[stat_id].loc_id != loc_id)
1044 continue;
1045
1046 const char *name
1047 = cs_map_name_to_id_reverse(_name_map, _sstats[stat_id].name_id);
1048
1049 double t_weight = -1;
1050 if (total_weight > 0 && _sstats[stat_id].intensive)
1051 t_weight = total_weight;
1052
1053 _log_array_info(" ",
1054 name,
1055 max_name_width,
1056 _sstats[stat_id].dim,
1057 n_g_elts,
1058 t_weight,
1059 vmin + stat_id,
1060 vmax + stat_id,
1061 vsum + stat_id,
1062 wsum + stat_id,
1063 &fpe_flag);
1064
1065 } /* End of loop on stats */
1066
1067 } /* End of loop on mesh locations */
1068
1069 sstat_cat_start = sstat_cat_end;
1070
1071 } /* End of loop on mesh categories */
1072
1073 /* Check NaN and exit */
1074 if (fpe_flag == 1)
1075 bft_error(__FILE__, __LINE__, 0,
1076 _("Invalid (not-a-number) values detected for a statistic."));
1077
1078 BFT_FREE(wsum);
1079 BFT_FREE(vsum);
1080 BFT_FREE(vmax);
1081 BFT_FREE(vmin);
1082
1083 cs_log_printf(CS_LOG_DEFAULT, "\n");
1084 }
1085
1086 /*----------------------------------------------------------------------------
1087 * Add or update clipping info for a given array
1088 *
1089 * parameters:
1090 * name_id Associated name id if not a field, -1 for a field
1091 * f_id associated field id, or -1
1092 * dim associated dimension
1093 * n_clip_min number of local clippings to minimum value
1094 * n_clip_max number of local clippings to maximum value
1095 * min_pre_clip minimum values prior to clipping
1096 * max_pre_clip maximum values prior to clipping
1097 *----------------------------------------------------------------------------*/
1098
1099 static void
_add_clipping(int name_id,int f_id,int dim,cs_lnum_t n_clip_min,cs_lnum_t n_clip_max,const cs_real_t min_pre_clip[],const cs_real_t max_pre_clip[],cs_lnum_t n_clip_min_comp[],cs_lnum_t n_clip_max_comp[])1100 _add_clipping(int name_id,
1101 int f_id,
1102 int dim,
1103 cs_lnum_t n_clip_min,
1104 cs_lnum_t n_clip_max,
1105 const cs_real_t min_pre_clip[],
1106 const cs_real_t max_pre_clip[],
1107 cs_lnum_t n_clip_min_comp[],
1108 cs_lnum_t n_clip_max_comp[])
1109 {
1110 bool need_sort = false;
1111
1112 int clip_id = _find_clip(f_id, name_id);
1113
1114 /* If not found, insert statistic */
1115
1116 if (clip_id < 0) {
1117
1118 _n_clips += 1;
1119 _clips_val_size += (dim == 1) ? 1 : dim + 1;
1120
1121 /* Reallocate key definitions if necessary */
1122
1123 if (_n_clips > _n_clips_max) {
1124 if (_n_clips_max == 0)
1125 _n_clips_max = 1;
1126 else
1127 _n_clips_max *= 2;
1128 BFT_REALLOC(_clips, _n_clips_max, cs_log_clip_t);
1129 }
1130
1131 if (_clips_val_size > _clips_val_size_max) {
1132 if (_clips_val_size_max == 0)
1133 _clips_val_size_max = dim;
1134 while (_clips_val_size > _clips_val_size_max)
1135 _clips_val_size_max *= 2;
1136 BFT_REALLOC(_clips_vmin, _clips_val_size_max, double);
1137 BFT_REALLOC(_clips_vmax, _clips_val_size_max, double);
1138 BFT_REALLOC(_clips_count, _clips_val_size_max*2, cs_gnum_t);
1139 }
1140
1141 need_sort = true; /* allow for binary search */
1142
1143 clip_id = _n_clips - 1;
1144 _clips[clip_id].f_id = f_id;
1145 _clips[clip_id].name_id = name_id;
1146 _clips[clip_id].dim = dim;
1147 _clips[clip_id].v_idx = (dim == 1) ? _clips_val_size - dim : _clips_val_size - dim - 1;
1148
1149 }
1150
1151 if (_clips[clip_id].dim != dim) {
1152 if (f_id > -1)
1153 bft_error(__FILE__, __LINE__, 0,
1154 "Clipping of field id %d previously defined in %s\n"
1155 "with dimension %d, redefined with dimension %d.",
1156 f_id, __func__,
1157 _clips[clip_id].dim, dim);
1158 else
1159 bft_error(__FILE__, __LINE__, 0,
1160 "Clipping of name %s previously defined in %s\n"
1161 "with dimension %d, redefined with dimension %d.",
1162 cs_map_name_to_id_reverse(_name_map, name_id),
1163 __func__,
1164 _clips[clip_id].dim, dim);
1165 }
1166
1167 /* Update clips */
1168
1169 _clips[clip_id].dim = dim;
1170
1171 int v_idx = _clips[clip_id].v_idx;
1172
1173 /* Prepare for future binary search */
1174
1175 if (need_sort)
1176 qsort(_clips, _n_clips, sizeof(cs_log_clip_t), &_compare_clips);
1177
1178 /* Update values */
1179 if (dim > 1) {
1180 _clips_count[(v_idx)*2] = n_clip_min;
1181 _clips_count[(v_idx)*2 + 1] = n_clip_max;
1182 _clips_vmin[v_idx] = min_pre_clip[0];
1183 _clips_vmax[v_idx] = max_pre_clip[0];
1184 for (int i = 0; i < dim; i++) {
1185 _clips_vmin[v_idx + i + 1] = min_pre_clip[i];
1186 _clips_vmax[v_idx + i + 1] = max_pre_clip[i];
1187 _clips_count[(v_idx + i + 1)*2] = n_clip_min_comp[i];
1188 _clips_count[(v_idx + i + 1)*2 + 1] = n_clip_max_comp[i];
1189 }
1190 }
1191 else {
1192 _clips_vmin[v_idx] = min_pre_clip[0];
1193 _clips_vmax[v_idx] = max_pre_clip[0];
1194 _clips_count[(v_idx)*2] = n_clip_min;
1195 _clips_count[(v_idx)*2 + 1] = n_clip_max;
1196 }
1197 }
1198
1199 /*----------------------------------------------------------------------------
1200 * Main logging output of additional clippings
1201 *----------------------------------------------------------------------------*/
1202
1203 static void
_log_clips(void)1204 _log_clips(void)
1205 {
1206 int clip_id;
1207 int type_idx[] = {0, 0, 0};
1208 double *vmin = NULL, *vmax = NULL;
1209 cs_gnum_t *vcount = NULL;
1210 size_t max_name_width = cs_log_strlen(_("field"));
1211 const int label_key_id = cs_field_key_id("label");
1212
1213 char tmp_s[5][64] = {"", "", "", "", ""};
1214
1215 const char *_cat_name[] = {N_("field"), N_("value")};
1216 const char *_cat_prefix[] = {"a ", "a "};
1217
1218 /* Allocate working arrays */
1219
1220 BFT_MALLOC(vmin, _clips_val_size, double);
1221 BFT_MALLOC(vmax, _clips_val_size, double);
1222 BFT_MALLOC(vcount, _clips_val_size*2, cs_gnum_t);
1223
1224 memcpy(vmin, _clips_vmin, _clips_val_size*sizeof(double));
1225 memcpy(vmax, _clips_vmax, _clips_val_size*sizeof(double));
1226 memcpy(vcount, _clips_count, _clips_val_size*sizeof(cs_gnum_t)*2);
1227
1228 /* Group MPI operations if required */
1229
1230 cs_parall_min(_clips_val_size, CS_DOUBLE, vmin);
1231 cs_parall_max(_clips_val_size, CS_DOUBLE, vmax);
1232 cs_parall_sum(_clips_val_size*2, CS_GNUM_TYPE, vcount);
1233
1234 /* Fist loop on clippings for counting */
1235
1236 for (clip_id = 0; clip_id < _n_clips; clip_id++) {
1237
1238 const char *name = NULL;
1239 int f_id = _clips[clip_id].f_id;
1240 int f_dim = 0;
1241 if (f_id > -1) {
1242 const cs_field_t *f = cs_field_by_id(f_id);
1243 name = cs_field_get_key_str(f, label_key_id);
1244 if (name == NULL)
1245 name = f->name;
1246 type_idx[1] = clip_id + 1;
1247 f_dim = f->dim;
1248 }
1249 else {
1250 name = cs_map_name_to_id_reverse(_name_map, _clips[clip_id].name_id);
1251 type_idx[2] = clip_id + 1;
1252 }
1253
1254 assert(name != NULL);
1255
1256 size_t l_name_width = cs_log_strlen(name);
1257 if (f_dim == 3)
1258 l_name_width += 3;
1259 else if (f_dim > 3)
1260 l_name_width += 4;
1261 max_name_width = CS_MAX(max_name_width, l_name_width);
1262
1263 }
1264
1265 if (type_idx[2] - type_idx[1] > 0) {
1266 size_t v_name_w = cs_log_strlen(_("value"));
1267 max_name_width = CS_MAX(max_name_width, v_name_w);
1268 }
1269
1270 max_name_width = CS_MIN(max_name_width, 63);
1271
1272 /* Loop on types */
1273
1274 for (int cat_id = 0; cat_id < 2; cat_id++) {
1275
1276 int start_id = type_idx[cat_id];
1277 int end_id = type_idx[cat_id+1];
1278
1279 if (end_id - start_id < 1)
1280 continue;
1281
1282 /* Print headers */
1283
1284 if (cat_id == 0)
1285 cs_log_printf(CS_LOG_DEFAULT,
1286 _("\n"
1287 " ** Clippings for computed fields\n"
1288 " -----------------------------\n"));
1289 else
1290 cs_log_printf(CS_LOG_DEFAULT,
1291 _("\n"
1292 " ** Clippings for auxiliary values\n"
1293 " ------------------------------\n"));
1294
1295 cs_log_strpad(tmp_s[0], _(_cat_name[cat_id]), max_name_width, 64);
1296 cs_log_strpadl(tmp_s[1], _("initial min"), 14, 64);
1297 cs_log_strpadl(tmp_s[2], _("initial max"), 14, 64);
1298 cs_log_strpadl(tmp_s[3], _("clips to min"), 12, 64);
1299 cs_log_strpadl(tmp_s[4], _("clips to max"), 12, 64);
1300 cs_log_printf(CS_LOG_DEFAULT,
1301 "\n %s %s %s %s %s\n",
1302 tmp_s[0], tmp_s[1], tmp_s[2], tmp_s[3], tmp_s[4]);
1303
1304 /* Underline */
1305
1306 for (size_t col = 0; col < 5; col++) {
1307 size_t i;
1308 size_t w0;
1309 if (col == 0)
1310 w0 = max_name_width;
1311 else if (col < 3)
1312 w0 = 14;
1313 else
1314 w0 = 12;
1315 for (i = 0; i < w0; i++)
1316 tmp_s[col][i] = '-';
1317 tmp_s[col][w0] = '\0';
1318 }
1319 cs_log_printf(CS_LOG_DEFAULT,
1320 "- %s %s %s %s %s\n",
1321 tmp_s[0], tmp_s[1], tmp_s[2], tmp_s[3], tmp_s[4]);
1322
1323 /* Second loop on clippings */
1324
1325 for (clip_id = start_id; clip_id < end_id; clip_id++) {
1326
1327 int v_idx = _clips[clip_id].v_idx;
1328
1329 const char *name = NULL;
1330 int f_id = _clips[clip_id].f_id;
1331 if (f_id > -1) {
1332 const cs_field_t *f = cs_field_by_id(f_id);
1333 name = cs_field_get_key_str(f, label_key_id);
1334 if (name == NULL)
1335 name = f->name;
1336 }
1337 else
1338 name = cs_map_name_to_id_reverse(_name_map, _clips[clip_id].name_id);
1339
1340 assert(name != NULL);
1341
1342 _log_clip_info(_cat_prefix[cat_id],
1343 name,
1344 max_name_width,
1345 _clips[clip_id].dim,
1346 vcount + v_idx*2,
1347 vcount + v_idx*2 + 1,
1348 vmin + v_idx,
1349 vmax + v_idx);
1350 }
1351
1352 }
1353
1354 BFT_FREE(vcount);
1355 BFT_FREE(vmax);
1356 BFT_FREE(vmin);
1357
1358 cs_log_printf(CS_LOG_DEFAULT, "\n");
1359 }
1360
1361 /*============================================================================
1362 * Fortran wrapper function definitions
1363 *============================================================================*/
1364
1365 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
1366
1367 /*============================================================================
1368 * Public function definitions
1369 *============================================================================*/
1370
1371 /*----------------------------------------------------------------------------*/
1372 /*!
1373 * \brief Free arrays possible used by logging of array statistics.
1374 */
1375 /*----------------------------------------------------------------------------*/
1376
1377 void
cs_log_iteration_destroy_all(void)1378 cs_log_iteration_destroy_all(void)
1379 {
1380 if (_category_map != NULL) {
1381 _sstats_val_size = 0;
1382 _sstats_val_size_max = 0;
1383 _n_sstats = 0;
1384 _n_sstats_max = 0;
1385 BFT_FREE(_sstats_vmin);
1386 BFT_FREE(_sstats_vmax);
1387 BFT_FREE(_sstats_vsum);
1388 BFT_FREE(_sstats_wsum);
1389 BFT_FREE(_sstats);
1390 cs_map_name_to_id_destroy(&_category_map);
1391 }
1392
1393 if (_n_clips_max > 0) {
1394 _clips_val_size = 0;
1395 _clips_val_size_max = 0;
1396 _n_clips = 0;
1397 _n_clips_max = 0;
1398 BFT_FREE(_clips_count);
1399 BFT_FREE(_clips_vmin);
1400 BFT_FREE(_clips_vmax);
1401 BFT_FREE(_clips);
1402 }
1403
1404 if (_name_map != NULL)
1405 cs_map_name_to_id_destroy(&_name_map);
1406
1407 if (_l2_residual_plot != NULL)
1408 cs_time_plot_finalize(&_l2_residual_plot);
1409 }
1410
1411 /*----------------------------------------------------------------------------*/
1412 /*!
1413 * \brief Log field and other array statistics for the current time step.
1414 */
1415 /*----------------------------------------------------------------------------*/
1416
1417 void
cs_log_iteration(void)1418 cs_log_iteration(void)
1419 {
1420 if (_n_clips > 0)
1421 _log_clips();
1422
1423 _log_fields();
1424
1425 if (_n_sstats > 0)
1426 _log_sstats();
1427
1428 cs_time_moment_log_iteration();
1429 cs_lagr_stat_log_iteration();
1430 cs_lagr_log_iteration();
1431
1432 cs_fan_log_iteration();
1433 cs_ctwr_log_balance();
1434 }
1435
1436 /*----------------------------------------------------------------------------*/
1437 /*!
1438 * \brief Add or update array not saved as permanent field to iteration log.
1439 *
1440 * \param[in] name array name
1441 * \param[in] category category name
1442 * \param[in] loc_id associated mesh location id
1443 * \param[in] is_intensive are the matching values intensive ?
1444 * \param[in] dim associated dimension (interleaved)
1445 * \param[in] val associated values
1446 */
1447 /*----------------------------------------------------------------------------*/
1448
1449 void
cs_log_iteration_add_array(const char * name,const char * category,const cs_mesh_location_type_t loc_id,bool is_intensive,int dim,const cs_real_t val[])1450 cs_log_iteration_add_array(const char *name,
1451 const char *category,
1452 const cs_mesh_location_type_t loc_id,
1453 bool is_intensive,
1454 int dim,
1455 const cs_real_t val[])
1456 {
1457 /* Initialize if necessary */
1458
1459 if (_name_map == NULL)
1460 _name_map = cs_map_name_to_id_create();
1461
1462 if (_category_map == NULL)
1463 _category_map = cs_map_name_to_id_create();
1464
1465 /* Find or insert entries in map */
1466
1467 bool need_sort = false;
1468
1469 int cat_id = cs_map_name_to_id(_category_map, category);
1470 int name_id = cs_map_name_to_id(_name_map, name);
1471
1472 int sstat_id = _find_sstats(cat_id, name_id);
1473
1474 /* If not found, insert statistic */
1475
1476 if (sstat_id < 0) {
1477
1478 int _dim = (dim == 3) ? 4 : dim;
1479
1480 _n_sstats += 1;
1481 _sstats_val_size += _dim;
1482
1483 /* Reallocate key definitions if necessary */
1484
1485 if (_n_sstats > _n_sstats_max) {
1486 if (_n_sstats_max == 0)
1487 _n_sstats_max = 1;
1488 else
1489 _n_sstats_max *= 2;
1490 BFT_REALLOC(_sstats, _n_sstats_max, cs_log_sstats_t);
1491 }
1492
1493 if (_sstats_val_size > _sstats_val_size_max) {
1494 if (_sstats_val_size_max == 0)
1495 _sstats_val_size_max = dim;
1496 while (_sstats_val_size > _sstats_val_size_max)
1497 _sstats_val_size_max *= 2;
1498 BFT_REALLOC(_sstats_vmin, _sstats_val_size_max, double);
1499 BFT_REALLOC(_sstats_vmax, _sstats_val_size_max, double);
1500 BFT_REALLOC(_sstats_vsum, _sstats_val_size_max, double);
1501 BFT_REALLOC(_sstats_wsum, _sstats_val_size_max, double);
1502 }
1503
1504 need_sort = true; /* allow for binary search */
1505
1506 sstat_id = _n_sstats - 1;
1507 _sstats[sstat_id].name_id = name_id;
1508 _sstats[sstat_id].cat_id = cat_id;
1509 _sstats[sstat_id].dim = dim;
1510 _sstats[sstat_id].v_idx = _sstats_val_size - _dim;
1511
1512 }
1513
1514 if (_sstats[sstat_id].dim != dim)
1515 bft_error(__FILE__, __LINE__, 0,
1516 "Array of name %s and category %s previously defined in %s\n"
1517 "with dimension %d, redefined with dimension %d.",
1518 cs_map_name_to_id_reverse(_name_map, name_id),
1519 cs_map_name_to_id_reverse(_category_map, cat_id),
1520 __func__,
1521 _sstats[sstat_id].dim, dim);
1522
1523 /* Update stats */
1524
1525 _sstats[sstat_id].loc_id = loc_id;
1526 _sstats[sstat_id].intensive = is_intensive;
1527 _sstats[sstat_id].dim = dim;
1528
1529 int v_idx = _sstats[sstat_id].v_idx;
1530
1531 /* Prepare for future binary search */
1532
1533 if (need_sort)
1534 qsort(_sstats, _n_sstats, sizeof(cs_log_sstats_t), &_compare_sstats);
1535
1536 /* Compute sstats */
1537
1538 bool have_weight = false;
1539 const cs_real_t *weight = NULL;
1540
1541 if (is_intensive != false) {
1542 cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
1543 switch(loc_id) {
1544 case CS_MESH_LOCATION_CELLS:
1545 weight = mq->cell_vol;
1546 have_weight = true;
1547 break;
1548 case CS_MESH_LOCATION_INTERIOR_FACES:
1549 weight = mq->i_face_surf;
1550 have_weight = true;
1551 break;
1552 case CS_MESH_LOCATION_BOUNDARY_FACES:
1553 weight = mq->b_face_surf;
1554 have_weight = true;
1555 break;
1556 case CS_MESH_LOCATION_VERTICES:
1557 have_weight = false;
1558 break;
1559 default:
1560 break;
1561 }
1562 }
1563
1564 const cs_lnum_t *n_elts = cs_mesh_location_get_n_elts(loc_id);
1565 const cs_lnum_t *elt_list = cs_mesh_location_get_elt_list(loc_id);
1566
1567 if (have_weight)
1568 cs_array_reduce_simple_stats_l_w(n_elts[0],
1569 dim,
1570 elt_list,
1571 elt_list,
1572 val,
1573 weight,
1574 _sstats_vmin + v_idx,
1575 _sstats_vmax + v_idx,
1576 _sstats_vsum + v_idx,
1577 _sstats_wsum + v_idx);
1578 else {
1579 cs_array_reduce_simple_stats_l(n_elts[0],
1580 dim,
1581 elt_list,
1582 val,
1583 _sstats_vmin + v_idx,
1584 _sstats_vmax + v_idx,
1585 _sstats_vsum + v_idx);
1586 for (int i = 0; i < dim; i++)
1587 _sstats_wsum[v_idx + i] = 0;
1588 }
1589 }
1590
1591 /*----------------------------------------------------------------------------*/
1592 /*!
1593 * \brief Add or update clipping info for a given array.
1594 *
1595 * \param[in] name field or array name
1596 * \param[in] dim associated dimension
1597 * \param[in] n_clip_min number of local clippings to minimum value
1598 * \param[in] n_clip_max number of local clippings to maximum value
1599 * \param[in] min_pre_clip minimum values prior to clipping
1600 * \param[in] max_pre_clip maximum values prior to clipping
1601 */
1602 /*----------------------------------------------------------------------------*/
1603
1604 void
cs_log_iteration_clipping(const char * name,int dim,cs_lnum_t n_clip_min,cs_lnum_t n_clip_max,const cs_real_t min_pre_clip[],const cs_real_t max_pre_clip[])1605 cs_log_iteration_clipping(const char *name,
1606 int dim,
1607 cs_lnum_t n_clip_min,
1608 cs_lnum_t n_clip_max,
1609 const cs_real_t min_pre_clip[],
1610 const cs_real_t max_pre_clip[])
1611 {
1612 /* Initialize if necessary */
1613
1614 if (_name_map == NULL)
1615 _name_map = cs_map_name_to_id_create();
1616
1617 int name_id = cs_map_name_to_id(_name_map, name);
1618
1619 _add_clipping(name_id, -1, dim,
1620 n_clip_min, n_clip_max,
1621 min_pre_clip, max_pre_clip, 0, 0);
1622 }
1623
1624 /*----------------------------------------------------------------------------*/
1625 /*!
1626 * \brief Add or update clipping info for a field
1627 *
1628 * \param[in] f_id associated field id
1629 * \param[in] n_clip_min number of local clippings to minimum value
1630 * \param[in] n_clip_max number of local clippings to maximum value
1631 * \param[in] min_pre_clip minimum values prior to clipping
1632 * \param[in] max_pre_clip maximum values prior to clipping
1633 * \param[in] n_clip_min_comp number of clip min by component
1634 * \param[in] n_clip_max_comp number of clip max by component
1635 */
1636 /*----------------------------------------------------------------------------*/
1637
1638 void
cs_log_iteration_clipping_field(int f_id,cs_lnum_t n_clip_min,cs_lnum_t n_clip_max,const cs_real_t min_pre_clip[],const cs_real_t max_pre_clip[],cs_lnum_t n_clip_min_comp[],cs_lnum_t n_clip_max_comp[])1639 cs_log_iteration_clipping_field(int f_id,
1640 cs_lnum_t n_clip_min,
1641 cs_lnum_t n_clip_max,
1642 const cs_real_t min_pre_clip[],
1643 const cs_real_t max_pre_clip[],
1644 cs_lnum_t n_clip_min_comp[],
1645 cs_lnum_t n_clip_max_comp[])
1646 {
1647 const cs_field_t *f = cs_field_by_id(f_id);
1648
1649 _add_clipping(-1, f_id, f->dim,
1650 n_clip_min, n_clip_max,
1651 min_pre_clip, max_pre_clip,n_clip_min_comp,n_clip_max_comp);
1652 }
1653
1654 /*----------------------------------------------------------------------------*/
1655 /*!
1656 * \brief Log L2 time residual for every variable fields.
1657 */
1658 /*----------------------------------------------------------------------------*/
1659
1660 void
cs_log_l2residual(void)1661 cs_log_l2residual(void)
1662 {
1663 if (cs_glob_rank_id > 0)
1664 return;
1665
1666 const cs_time_step_t *ts = cs_glob_time_step;
1667 const int n_fields = cs_field_n_fields();
1668
1669 /* write header */
1670
1671 if (_l2_residual_plot == NULL) {
1672
1673 int _plot_buffer_steps = -1;
1674 double _plot_flush_wtime = 3600;
1675 cs_time_plot_format_t _plot_format = CS_TIME_PLOT_CSV;
1676 bool use_iteration = (ts->is_local) ? true : false;
1677
1678 const char **labels;
1679 BFT_MALLOC(labels, n_fields + 1, const char *);
1680
1681 int n_variables = 0;
1682 for (int f_id = 0 ; f_id < n_fields ; f_id++) {
1683 const cs_field_t *f = cs_field_by_id(f_id);
1684 if (f->type & CS_FIELD_VARIABLE) {
1685 labels[n_variables] = f->name;
1686 n_variables++;
1687 }
1688 }
1689
1690 _l2_residual_plot = cs_time_plot_init_probe("residuals",
1691 "",
1692 _plot_format,
1693 use_iteration,
1694 _plot_flush_wtime,
1695 _plot_buffer_steps,
1696 n_variables,
1697 NULL,
1698 NULL,
1699 labels);
1700
1701 BFT_FREE(labels);
1702 }
1703
1704 {
1705 cs_real_t *vals;
1706 BFT_MALLOC(vals, n_fields, cs_real_t);
1707
1708 int si_k_id = cs_field_key_id("solving_info");
1709
1710 int n_variables = 0;
1711 for (int f_id = 0 ; f_id < n_fields ; f_id++) {
1712 const cs_field_t *f = cs_field_by_id(f_id);
1713 if (f->type & CS_FIELD_VARIABLE) {
1714 const cs_solving_info_t *sinfo
1715 = cs_field_get_key_struct_const_ptr(f, si_k_id);
1716 vals[n_variables] = sinfo->l2residual;
1717 n_variables += 1;
1718 }
1719 }
1720
1721 cs_time_plot_vals_write(_l2_residual_plot,
1722 ts->nt_cur,
1723 ts->t_cur,
1724 n_variables,
1725 vals);
1726
1727 BFT_FREE(vals);
1728 }
1729 }
1730
1731 /*----------------------------------------------------------------------------*/
1732
1733 END_C_DECLS
1734