1 /*============================================================================
2 * Checkpoint/restart handling for default application.
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
31 /*----------------------------------------------------------------------------
32 * Standard C library headers
33 *----------------------------------------------------------------------------*/
34
35 #include <assert.h>
36 #include <ctype.h>
37 #include <math.h>
38 #include <stdio.h>
39 #include <stdlib.h>
40 #include <string.h>
41
42 /*----------------------------------------------------------------------------
43 * Local headers
44 *----------------------------------------------------------------------------*/
45
46 #include "bft_mem.h"
47 #include "bft_error.h"
48 #include "bft_printf.h"
49
50 #include "cs_array.h"
51 #include "cs_base.h"
52 #include "cs_field.h"
53 #include "cs_field_pointer.h"
54 #include "cs_field_operator.h"
55 #include "cs_gradient.h"
56 #include "cs_halo.h"
57 #include "cs_halo_perio.h"
58 #include "cs_les_inflow.h"
59 #include "cs_log.h"
60 #include "cs_map.h"
61 #include "cs_math.h"
62 #include "cs_mesh.h"
63 #include "cs_parall.h"
64 #include "cs_mesh_location.h"
65 #include "cs_random.h"
66 #include "cs_time_step.h"
67 #include "cs_turbulence_model.h"
68 #include "cs_physical_constants.h"
69 #include "cs_volume_zone.h"
70
71 /*----------------------------------------------------------------------------
72 * Header for the current file
73 *----------------------------------------------------------------------------*/
74
75 #include "cs_restart_default.h"
76
77 /*----------------------------------------------------------------------------*/
78
79 BEGIN_C_DECLS
80
81 /*=============================================================================
82 * Additional doxygen documentation
83 *============================================================================*/
84
85 /*!
86 \file cs_restart_default.c
87 Checkpoint/restart handling for default application.
88 */
89
90 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
91
92 /*=============================================================================
93 * Macro definitions
94 *============================================================================*/
95
96 /*============================================================================
97 * Prototypes for functions intended for use only by Fortran wrappers.
98 * (descriptions follow, with function bodies).
99 *============================================================================*/
100
101 /*============================================================================
102 * Type definitions
103 *============================================================================*/
104
105 /*============================================================================
106 * Static global variables
107 *============================================================================*/
108
109 const char *_coeff_name[] = {"bc_coeffs::a", "bc_coeffs::b",
110 "bc_coeffs::af", "bc_coeffs::bf",
111 "bc_coeffs::ad", "bc_coeffs::bd",
112 "bc_coeffs::ac", "bc_coeffs::bc"};
113
114 /*============================================================================
115 * Private function definitions
116 *============================================================================*/
117
118 /*----------------------------------------------------------------------------
119 * Cancel array values in solid zone
120 *
121 * parameters:
122 * stride <-- array stride
123 * a <-> array
124 *----------------------------------------------------------------------------*/
125
126 /*----------------------------------------------------------------------------
127 * Define internal coupling through the GUI.
128 *----------------------------------------------------------------------------*/
129
130 static void
_cancel_in_solid_zones(cs_lnum_t stride,cs_real_t * a)131 _cancel_in_solid_zones(cs_lnum_t stride,
132 cs_real_t *a)
133 {
134 int n_zones = cs_volume_zone_n_zones();
135
136 for (int i = 0; i < n_zones; i++) {
137 const cs_zone_t *z = cs_volume_zone_by_id(i);
138 if (z->type & CS_VOLUME_ZONE_SOLID) {
139 if (stride == 1) {
140 for (cs_lnum_t j = 0; j < z->n_elts; z++)
141 a[z->elt_ids[j]] = 0;
142 }
143 else {
144 for (cs_lnum_t j = 0; j < z->n_elts; z++) {
145 cs_real_t *_a = a + (stride * z->elt_ids[j]);
146 for (cs_lnum_t k = 0; k < stride; k++)
147 _a[k] = 0;
148 }
149 }
150 }
151 }
152 }
153 /*----------------------------------------------------------------------------
154 * Read and rebuild partial field metadata from legacy checkpoint.
155 *
156 * Note that when reading legacy files (Code_Saturne version 3.3 and below),
157 * the old id will actually be the old scalar id (-1 for others).
158 *
159 * parameters:
160 * r <-- associated restart file pointer
161 *----------------------------------------------------------------------------*/
162
163 static void
_read_legacy_field_info(cs_restart_t * r)164 _read_legacy_field_info(cs_restart_t *r)
165 {
166 int retcode;
167
168 int n_fields = cs_field_n_fields();
169
170 /* Initialization */
171
172 int kold = cs_field_key_id_try("old_scalar_num");
173
174 /* Now read headers */
175
176 cs_lnum_t n_old[4] = {0, 0, 0, 0}, n_cur[4] = {0, 0, 0, 0};
177
178 const char *sec_id[] = {"nombre_variables",
179 "nombre_scalaires",
180 "nombre_scalaires_us",
181 "nombre_scalaires_pp"};
182
183 for (int i = 0; i < 4; i++) {
184 retcode = cs_restart_read_section(r,
185 sec_id[i],
186 CS_MESH_LOCATION_NONE,
187 1,
188 CS_TYPE_int,
189 n_old + i);
190 if (retcode != CS_RESTART_SUCCESS)
191 bft_error
192 (__FILE__, __LINE__, 0,
193 _("Error reading variable information in restart file \"%s\"."),
194 cs_restart_get_name(r));
195 }
196
197 const int kv = cs_field_key_id_try("variable_id");
198 const int ks = cs_field_key_id_try("scalar_id");
199
200 /* Count variables and user and non-user scalars */
201
202 for (int f_id = 0; f_id < n_fields; f_id++) {
203 const cs_field_t *f = cs_field_by_id(f_id);
204 int v_num = -1;
205 if (kv > -1)
206 v_num = cs_field_get_key_int(f, kv);
207 if (v_num > 0) {
208 int s_num = -1;
209 n_cur[0] += 1;
210 if (ks > -1)
211 s_num = cs_field_get_key_int(f, ks);
212 if (s_num > -1) {
213 n_cur[1] += 1;
214 if (f->type & CS_FIELD_USER)
215 n_cur[2] += 1;
216 else
217 n_cur[3] += 1;
218 }
219 }
220 }
221
222 /* Possible shift in old ids if temperature has been moved from
223 user to model scalar */
224
225 int us_shift = 0, pp_shift = 0;
226
227 /* Special case if temperature has been moved from
228 user to model scalar */
229
230 if ( n_cur[1] == n_old[1] && n_cur[2] == n_old[2] -1
231 && n_cur[3] == 1 && n_old[3] == 0) {
232
233 if (CS_F_(t) != NULL || CS_F_(h) != NULL) {
234 us_shift = -1;
235 pp_shift = n_cur[2];
236 }
237
238 }
239
240 /* Warn in case of change */
241
242 if ( n_cur[0] != n_old[0] || n_cur[1] != n_old[1]
243 || n_cur[2] != n_old[2] || n_cur[3] != n_old[3]) {
244
245 /* Special case if temperature has been moved from
246 user to model scalar */
247
248 if (n_cur[0] == n_old[0] && n_cur[1] == n_old[1] && us_shift == -1)
249 bft_printf
250 (_("\nRemark: the thermal scalar was treated as a user scalar\n"
251 " in the restart file, and is moved to a model scalar\n"
252 " in the current computation.\n"));
253
254 else {
255 bft_printf
256 (_("\n"
257 " Warning: the number of variables or scalars has been changed\n"
258 " relative to the restart file.\n\n"
259 " currently %d variables, of which %d scalars\n"
260 " previously %d variables, of which %d scalars\n\n"
261 " The computation continues, with a partial restart.\n"),
262 (int)(n_cur[0]), (int)(n_cur[1]), (int)(n_old[0]), (int)(n_old[1]));
263
264 }
265
266 }
267
268 /* Now check (and update if necessary) old scalar id */
269
270 for (int f_id = 0; f_id < n_fields; f_id++) {
271 cs_field_t *f = cs_field_by_id(f_id);
272 int old_scal_num = -1;
273 int s_num = -1;
274 if (ks > -1)
275 s_num = cs_field_get_key_int(f, ks);
276 if (s_num > -1) {
277 old_scal_num = -1;
278 if (kold > -1)
279 old_scal_num = cs_field_get_key_int(f, kold);
280 if (old_scal_num < 0) {
281 if (f->type & CS_FIELD_USER)
282 old_scal_num = s_num + us_shift;
283 else
284 old_scal_num = s_num + pp_shift;
285 if (old_scal_num > n_old[1])
286 old_scal_num = -1;
287 }
288 else {
289 if (old_scal_num > n_old[1])
290 bft_error
291 (__FILE__, __LINE__, 0,
292 _("Field \"%s\" has user-defined key \"old_scalar_num\" value %d,\n"
293 "but the number of available scalars in restart is %d."),
294 f->name, old_scal_num, (int)(n_old[1]));
295 }
296 if (kold < 0)
297 kold = cs_field_define_key_int("old_scalar_num",
298 -1,
299 CS_FIELD_VARIABLE);
300 cs_field_set_key_int(f, kold, old_scal_num);
301 }
302 }
303
304 }
305
306 /*----------------------------------------------------------------------------
307 * Synchronize cell-based field values.
308 *
309 * parameters:
310 * f <-> field whose values should be synchronized
311 * t_id <-- time id (0 for current, 1 for previous, ...)
312 *----------------------------------------------------------------------------*/
313
314 static void
_sync_field_vals(cs_field_t * f,int t_id)315 _sync_field_vals(cs_field_t *f,
316 int t_id)
317 {
318 const cs_mesh_t *m = cs_glob_mesh;
319
320 if (m->halo != NULL) {
321
322 cs_halo_type_t halo_type = CS_HALO_EXTENDED;
323 cs_real_t *v = f->vals[t_id];
324
325 cs_halo_sync_var_strided(m->halo, halo_type, v, f->dim);
326
327 if (m->n_init_perio > 0) {
328 if (f->dim == 3)
329 cs_halo_perio_sync_var_vect(m->halo, halo_type, v, 3);
330 else if (f->dim == 6)
331 cs_halo_perio_sync_var_sym_tens(m->halo, halo_type, v);
332 else if (f->dim == 9)
333 cs_halo_perio_sync_var_tens(m->halo, halo_type, v);
334 }
335
336 }
337 }
338
339 /*----------------------------------------------------------------------------
340 * Read field values from checkpoint.
341 *
342 * Values are found using the default rules based on the field's name
343 * postfixed by ::vals::%t_id, or its name itself.
344 *
345 * parameters:
346 * r <-- associated restart file pointer
347 * restart_name <-- base name with which read is attempted
348 * t_id <-- time id (0 for current, 1 for previous, ...)
349 * f <-> field whose values should be read
350 *
351 * returns:
352 * CS_RESTART_SUCCESS in case of success, CS_RESTART_ERR_... otherwise
353 *----------------------------------------------------------------------------*/
354
355 static int
_read_field_vals(cs_restart_t * r,const char * r_name,int t_id,cs_field_t * f)356 _read_field_vals(cs_restart_t *r,
357 const char *r_name,
358 int t_id,
359 cs_field_t *f)
360 {
361 int retcode = CS_RESTART_SUCCESS;
362
363 char _sec_name[128];
364 char *sec_name = _sec_name;
365
366 if (strlen(r_name) > 96)
367 BFT_MALLOC(sec_name, strlen(r_name) + 64, char); /* wide margin */
368
369 /* Check for data; data will be read later, so that compatibility
370 checks may be done first; we really try reading the data only
371 at the end, so if it is not found, a warning will be logged only
372 once (and not once per test), and preferentially use the
373 base (non-compatibility) name. */
374
375 snprintf(sec_name, 127, "%s::vals::%d", r_name, t_id);
376 sec_name[127] = '\0';
377
378 retcode = cs_restart_check_section(r,
379 sec_name,
380 f->location_id,
381 f->dim,
382 CS_TYPE_cs_real_t);
383
384 /* Otherwise, try reading with basic (restart) name only if requested */
385
386 if ( (retcode == CS_RESTART_ERR_EXISTS || retcode == CS_RESTART_ERR_N_VALS)
387 && r_name != f->name) {
388 snprintf(sec_name, 127, "%s", r_name);
389 sec_name[127] = '\0';
390 retcode = cs_restart_check_section(r,
391 sec_name,
392 f->location_id,
393 f->dim,
394 CS_TYPE_cs_real_t);
395 }
396
397 /* Read if available */
398
399 if (retcode == CS_RESTART_SUCCESS)
400 retcode = cs_restart_read_section(r,
401 sec_name,
402 f->location_id,
403 f->dim,
404 CS_TYPE_cs_real_t,
405 f->vals[t_id]);
406
407 /* Try to read anyways (with base name) to log warning */
408
409 else {
410 snprintf(sec_name, 127, "%s::vals::%d", r_name, t_id);
411 sec_name[127] = '\0';
412 retcode = cs_restart_read_section(r,
413 sec_name,
414 f->location_id,
415 f->dim,
416 CS_TYPE_cs_real_t,
417 f->vals[t_id]);
418 }
419
420 if (sec_name != _sec_name)
421 BFT_FREE(sec_name);
422
423 if ( retcode == CS_RESTART_SUCCESS
424 && f->location_id == CS_MESH_LOCATION_CELLS)
425 _sync_field_vals(f, t_id);
426
427 return retcode;
428 }
429
430 /*----------------------------------------------------------------------------
431 * Read Rij field using interleaved or non-interleaved (legacy) sections.
432 *
433 * parameters:
434 * r <-- associated restart file pointer
435 * location_id <-- mesh locartion id
436 * t_id <-- time id (0 for current, 1 for previous, ...)
437 * rij <-> Rij values for chosen time step
438 *
439 * returns:
440 * CS_RESTART_SUCCESS in case of success, CS_RESTART_ERR_... otherwise
441 *----------------------------------------------------------------------------*/
442
443 static int
_read_rij(cs_restart_t * r,int location_id,int t_id,cs_real_t rij[][6])444 _read_rij(cs_restart_t *r,
445 int location_id,
446 int t_id,
447 cs_real_t rij[][6])
448 {
449 int retcode = cs_restart_read_real_6_t_compat(r,
450 "rij::vals::0",
451 "r11::vals::0",
452 "r22::vals::0",
453 "r33::vals::0",
454 "r12::vals::0",
455 "r23::vals::0",
456 "r13::vals::0",
457 location_id,
458 rij);
459
460 if (retcode == CS_RESTART_ERR_EXISTS && t_id == 0) {
461
462 const cs_lnum_t n_cells = cs_glob_mesh->n_cells;
463 const char *old_names[] = {"R11", "R22", "R33", "R12", "R23", "R13"};
464
465 cs_real_t *v_tmp;
466 BFT_MALLOC(v_tmp, n_cells, cs_real_t);
467
468 for (cs_lnum_t j = 0; j < 6; j++) {
469 char old_sec_name[128];
470 snprintf(old_sec_name, 127, "%s_ce_phase01", old_names[j]);
471 old_sec_name[127] = '\0';
472
473 retcode = cs_restart_check_section(r,
474 old_sec_name,
475 location_id,
476 1,
477 CS_TYPE_cs_real_t);
478 if (retcode != CS_RESTART_SUCCESS)
479 break;
480
481 retcode = cs_restart_read_section(r,
482 old_sec_name,
483 location_id,
484 1,
485 CS_TYPE_cs_real_t,
486 v_tmp);
487
488 if (retcode != CS_RESTART_SUCCESS)
489 break;
490
491 for (cs_lnum_t i = 0; i < n_cells; i++)
492 rij[i][j] = v_tmp[i];
493 }
494
495 BFT_FREE(v_tmp);
496 }
497
498 return retcode;
499 }
500
501 /*----------------------------------------------------------------------------
502 * Read field values from legacy checkpoint.
503 *
504 * Values are found using older names for compatibility with older files.
505 * For cell-based fields, the old name base is appended automatically with
506 * "_ce_phase01", except for scalars, where the name uses a different scheme,
507 * based on "scalaire_ce_%04" % s_num.
508 *
509 * parameters:
510 * r <-- associated restart file pointer
511 * r_name <-- base name with which read is attempted
512 * t_id <-- time id (0 for current, 1 for previous, ...)
513 * f <-> file whose values should be read
514 *
515 * returns:
516 * CS_RESTART_SUCCESS in case of success, CS_RESTART_ERR_... otherwise
517 *----------------------------------------------------------------------------*/
518
519 static int
_read_field_vals_legacy(cs_restart_t * r,const char * r_name,int t_id,cs_field_t * f)520 _read_field_vals_legacy(cs_restart_t *r,
521 const char *r_name,
522 int t_id,
523 cs_field_t *f)
524 {
525 char sec_name[128] = "";
526 char old_name_x[96] = "", old_name_y[96] = "", old_name_z[96] = "";
527
528 int retcode = CS_RESTART_SUCCESS;
529
530 /* Check for renaming */
531
532 char old_name[96] = "";
533 int ks = cs_field_key_id_try("scalar_id");
534 int scalar_id = cs_field_get_key_int(f, ks);
535
536 /* Special case for scalars */
537
538 if (scalar_id > -1) {
539 if (r_name != f->name) {
540 const char *name = r_name;
541 while (*name != '\0' && !isdigit(*name))
542 name++;
543 scalar_id = atoi(name) - 1;
544 }
545 if (scalar_id > -1)
546 snprintf(old_name, 96, "%04d", scalar_id+1);
547 else
548 snprintf(old_name, 96, "%s", r_name);
549 }
550
551 /* Other fields may need specific renaming
552 (old_name for partial section name, sec_name for direct section name) */
553
554 else if (r_name == f->name) {
555 snprintf(old_name, sizeof(old_name), "%s", f->name);
556 if (f == CS_F_(vel)) {
557 if (t_id == 0)
558 strncpy(old_name, "vitesse", sizeof(old_name));
559 else if (t_id == 1)
560 strncpy(sec_name, "velocity_prev", sizeof(old_name));
561 }
562 else if (f == CS_F_(p))
563 strncpy(old_name, "pression", sizeof(old_name));
564 else if (f == CS_F_(rij))
565 strncpy(old_name, "Rij", sizeof(old_name));
566 else if (f == CS_F_(eps))
567 strncpy(old_name, "eps", sizeof(old_name));
568 else if (f == CS_F_(f_bar))
569 strncpy(old_name, "fb", sizeof(old_name));
570 else if (f == CS_F_(alp_bl)) {
571 /* Special case: "al" also possible here, depending on turbulence model;
572 check for either, with one test for main restart, the other for
573 the auxilairy restart */
574 int sec_code;
575 strncpy(old_name, "alp", sizeof(old_name));
576 sec_code = cs_restart_check_section(r, "al_ce_phase01",
577 1, 1, CS_TYPE_cs_real_t);
578 if (sec_code == CS_RESTART_SUCCESS)
579 strncpy(old_name, "al", sizeof(old_name));
580 else
581 sec_code = cs_restart_check_section(r, "fm_al_phase01",
582 0, 1, CS_TYPE_int);
583 if (sec_code == CS_RESTART_SUCCESS)
584 strncpy(old_name, "al", sizeof(old_name));
585 }
586 else if (f == CS_F_(nusa))
587 strncpy(old_name, "nusa", sizeof(old_name));
588 else if (f == CS_F_(mesh_u))
589 strncpy(old_name, "vit_maillage", sizeof(old_name));
590 else if (f == CS_F_(rho)) {
591 if (t_id == 0)
592 strncpy(old_name, "rho", sizeof(old_name));
593 else if (t_id == 1)
594 strncpy(old_name, "rho_old", sizeof(old_name));
595 }
596 else if (f == CS_F_(rho_b))
597 strncpy(sec_name, "rho_fb_phase01", 96);
598
599 else if (f == CS_F_(cp))
600 strncpy(old_name, "cp", sizeof(old_name));
601
602 else if (f == CS_F_(mu))
603 strncpy(old_name, "viscl", sizeof(old_name));
604 else if (f == CS_F_(mu_t))
605 strncpy(old_name, "visct", sizeof(old_name));
606
607 else if (f == CS_F_(t_b))
608 strncpy(old_name, "tparoi_fb", sizeof(old_name));
609 else if (f == CS_F_(qinci))
610 strncpy(old_name, "qincid_fb", sizeof(old_name));
611 else if (f == CS_F_(hconv))
612 strncpy(old_name, "hfconv_fb", sizeof(old_name));
613 else if (f == CS_F_(fconv))
614 strncpy(old_name, "flconv_fb", sizeof(old_name));
615
616 else if (strcmp(f->name, "dt") == 0)
617 strncpy(sec_name, "dt_variable_espace_ce", 127);
618
619 else if (strcmp(f->name, "dt") == 0)
620 strncpy(sec_name, "dt_variable_espace_ce", 127);
621
622 else if ( f->location_id == CS_MESH_LOCATION_VERTICES
623 && strcmp(f->name, "mesh_displacement") == 0
624 && t_id == 0)
625 strncpy(sec_name, "vertex_displacement", 127);
626
627 else if (strcmp(f->name, "void_fraction") == 0)
628 strncpy(sec_name, "taux_vide_ce", 127);
629
630 else if (strcmp(f->name, "rad_st") == 0)
631 strncpy(sec_name, "rayexp_ce", 127);
632 else if (strcmp(f->name, "rad_st_implicit") == 0)
633 strncpy(sec_name, "rayimp_ce", 127);
634 else if (f == CS_F_(rad_lumin))
635 strncpy(sec_name, "luminance", 127);
636
637 else if (strcmp(f->name, "joule_power") == 0)
638 strncpy(sec_name, "tsource_sc_ce_joule", 127);
639 else if (strcmp(f->name, "laplace_force") == 0)
640 strncpy(old_name, "laplace_force", sizeof(old_name));
641 }
642
643 if (sec_name[0] == '\0') {
644 if (scalar_id > -1)
645 snprintf(sec_name, 127, "scalaire_ce_%04d", scalar_id);
646 else if (f->location_id == CS_MESH_LOCATION_CELLS)
647 snprintf(sec_name, 127, "%s_ce_phase01", old_name);
648 else
649 snprintf(sec_name, 127, "%s", old_name);
650 }
651
652 sec_name[127] = '\0';
653
654 retcode = cs_restart_check_section(r,
655 sec_name,
656 f->location_id,
657 f->dim,
658 CS_TYPE_cs_real_t);
659
660 if (retcode == CS_RESTART_SUCCESS)
661 retcode = cs_restart_read_section(r,
662 sec_name,
663 f->location_id,
664 f->dim,
665 CS_TYPE_cs_real_t,
666 f->vals[t_id]);
667
668 /* Last chance for 3D fields */
669
670 else if (f->dim == 3 && retcode == CS_RESTART_ERR_EXISTS) {
671
672 if (strcmp(old_name, "vit_maillage") == 0) {
673 snprintf(old_name_x, 96, "%s_u_ce", old_name);
674 snprintf(old_name_y, 96, "%s_v_ce", old_name);
675 snprintf(old_name_z, 96, "%s_w_ce", old_name);
676 }
677 else if (strcmp(old_name, "laplace_force") == 0) {
678 snprintf(old_name_x, 96, "%s_1", old_name);
679 snprintf(old_name_y, 96, "%s_2", old_name);
680 snprintf(old_name_z, 96, "%s_2", old_name);
681 }
682 else {
683 snprintf(old_name_x, 96, "%s_u_ce_phase01", old_name);
684 snprintf(old_name_y, 96, "%s_v_ce_phase01", old_name);
685 snprintf(old_name_z, 96, "%s_w_ce_phase01", old_name);
686 }
687
688 old_name_x[95] = '\0';
689 old_name_y[95] = '\0';
690 old_name_z[95] = '\0';
691
692 retcode = cs_restart_check_section(r,
693 old_name_x,
694 f->location_id,
695 1,
696 CS_TYPE_cs_real_t);
697
698 if (retcode == CS_RESTART_SUCCESS)
699 retcode = cs_restart_read_real_3_t_compat(r,
700 sec_name,
701 old_name_x,
702 old_name_y,
703 old_name_z,
704 f->location_id,
705 (cs_real_3_t *)(f->vals[t_id]));
706 }
707 else if (f->dim == 6 && retcode == CS_RESTART_ERR_EXISTS) {
708
709 if (strcmp(old_name, "Rij") == 0) {
710 retcode = cs_restart_check_section(r,
711 "r11::vals::0",
712 f->location_id,
713 1,
714 CS_TYPE_cs_real_t);
715
716 if (retcode == CS_RESTART_SUCCESS)
717 _read_rij(r, f->location_id, t_id, (cs_real_6_t *)(f->vals[t_id]));
718 }
719 }
720
721 if ( retcode == CS_RESTART_SUCCESS
722 && f->location_id == CS_MESH_LOCATION_CELLS)
723 _sync_field_vals(f, t_id);
724
725 return retcode;
726 }
727
728 /*----------------------------------------------------------------------------
729 * Determine mass flux number of a given variable for legacy restart file
730 *
731 * parameters:
732 * r <-> associated restart file pointer
733 * f <-- associated field pointer
734 * scalar_num <-- associated scalar number, or -1
735 * t_id <-- associated time id (0: current, 1: previous)
736 *
737 * returns:
738 * number of matching mass flux in restart file (-1 if read failed)
739 *----------------------------------------------------------------------------*/
740
741 static int
_legacy_mass_flux_num(cs_restart_t * r,const cs_field_t * f,int scalar_num,int t_id)742 _legacy_mass_flux_num(cs_restart_t *r,
743 const cs_field_t *f,
744 int scalar_num,
745 int t_id)
746 {
747 int retval = 1;
748
749 /* As of Code_Saturne 3.3, only scalars may have a different mass flux
750 from the "main" mass flux (in the case of scalars with drift), so for
751 all others, reading the associated mass flux name is of no real use. */
752
753 char sec_name[128] = "";
754
755 const char *prefix[2] = {"fm_", "fm_a_"};
756 if (scalar_num > 0)
757 snprintf(sec_name, 127, "%sscalaire%04d", prefix[t_id], scalar_num);
758 else if (strcmp(f->name, "void_fraction") == 0)
759 snprintf(sec_name, 127, "%staux_vide", prefix[t_id]);
760
761 /* Read from restart */
762
763 if (sec_name[0] != '\0') {
764 cs_lnum_t buf[1];
765 sec_name[127] = '\0';
766 int retcode = cs_restart_read_section(r,
767 sec_name,
768 CS_MESH_LOCATION_NONE,
769 1,
770 CS_TYPE_int,
771 buf);
772 if (retcode == CS_RESTART_SUCCESS)
773 retval = buf[0];
774 else
775 retval = -1;
776 }
777
778 return retval;
779 }
780
781 /*----------------------------------------------------------------------------
782 * Read fields depending on others from checkpoint.
783 *
784 * This function handles legacy files (Code_Saturne version 3.3 and below).
785 *
786 * parameters:
787 * r <-> associated restart file pointer
788 * key <-> key for field association
789 * read_flag <-> flag to track fields read, or NULL;
790 * set to 1 for fields read, -1 for fields
791 * failed to read (size: n_fields)
792 *
793 * returns:
794 * number of fields read
795 *----------------------------------------------------------------------------*/
796
797 static int
_read_linked_fields_legacy(cs_restart_t * r,const char * key,int read_flag[])798 _read_linked_fields_legacy(cs_restart_t *r,
799 const char *key,
800 int read_flag[])
801 {
802 int retcode;
803
804 int retcount = 0;
805
806 /* Initialization */
807
808 int category = 0;
809
810 const int n_fields = cs_field_n_fields();
811
812 const int key_id = cs_field_key_id(key);
813 const int key_flag = cs_field_key_flag(key_id);
814
815 const int kold = cs_field_key_id_try("old_scalar_num");
816 const int ks = cs_field_key_id_try("scalar_id");
817
818 /* Determine field type (out of possibilities in legacy files) */
819
820 if (strcmp(key, "inner_mass_flux_id") == 0)
821 category = 1;
822 else if (strcmp(key, "boundary_mass_flux_id") == 0)
823 category = 2;
824 else if (strcmp(key, "diffusivity_id") == 0)
825 category = 3;
826
827 for (int f_id = 0; f_id < n_fields; f_id++) {
828
829 const cs_field_t *f = cs_field_by_id(f_id);
830
831 if (key_flag == -1 || !(f->type & key_flag))
832 continue;
833
834 const int lnk_f_id = cs_field_get_key_int(f, key_id);
835 int s_num = -1;
836
837 if (lnk_f_id > -1) {
838
839 cs_field_t *f_lnk = cs_field_by_id(lnk_f_id);
840
841 if (read_flag[lnk_f_id] != 0)
842 continue;
843
844 read_flag[lnk_f_id] = -1;
845
846 /* check for (possibly renumbered) scalar */
847 if (f->type & CS_FIELD_VARIABLE) {
848 if (kold > -1)
849 s_num = cs_field_get_key_int(f, kold);
850 if (s_num < 0 && ks > -1)
851 s_num = cs_field_get_key_int(f, ks);
852 }
853
854 for (int t_id = 0; t_id < 2; t_id++) {
855
856 if (t_id >= f_lnk->n_time_vals)
857 break;
858
859 /* Build field name to read based on category */
860
861 char sec_name[128];
862 if (category == 1) {
863 int mf_num = _legacy_mass_flux_num(r, f, s_num, t_id);
864 if (t_id == 0)
865 snprintf(sec_name, 127, "flux_masse_fi_%04d", mf_num);
866 else
867 snprintf(sec_name, 127, "flux_masse_a_fi_%04d", mf_num);
868 }
869 else if (category == 2) {
870 int mf_num = _legacy_mass_flux_num(r, f, s_num, t_id);
871 if (t_id == 0)
872 snprintf(sec_name, 127, "flux_masse_fb_%04d", mf_num);
873 else
874 snprintf(sec_name, 127, "flux_masse_a_fb_%04d", mf_num);
875 }
876 else if (category == 3)
877 snprintf(sec_name, 127, "visls_ce_scalaire%04d", s_num);
878
879 /* Now we know which field name to read */
880
881 retcode = cs_restart_check_section(r,
882 sec_name,
883 f->location_id,
884 f->dim,
885 CS_TYPE_cs_real_t);
886
887 if (retcode == CS_RESTART_SUCCESS)
888 retcode = cs_restart_read_section(r,
889 sec_name,
890 f->location_id,
891 f->dim,
892 CS_TYPE_cs_real_t,
893 f->vals[t_id]);
894
895 if (retcode == CS_RESTART_SUCCESS) {
896 if (t_id == 0)
897 read_flag[lnk_f_id] = 1;
898 else
899 read_flag[lnk_f_id] += 2;
900 retcount += 1;
901 }
902
903 } /* t_id */
904
905 }
906
907 }
908
909 return retcount;
910 }
911
912 /*----------------------------------------------------------------------------
913 * Compare old and new values for a given field key
914 *
915 * parameters:
916 * r <-- associated restart file pointer
917 * key <-- associated model key
918 * old_field_map <-- name to id map of fields in restart file
919 *
920 * returns:
921 * number of values changed, or -1 if information was not found
922 *----------------------------------------------------------------------------*/
923
924 static int
_check_field_model(cs_restart_t * r,const char * key,const cs_map_name_to_id_t * old_field_map)925 _check_field_model(cs_restart_t *r,
926 const char *key,
927 const cs_map_name_to_id_t *old_field_map)
928 {
929 int retcode;
930
931 const int n_fields = cs_field_n_fields();
932 const int n_o_fields = cs_map_name_to_id_size(old_field_map);
933
934 const int key_id = cs_field_key_id(key);
935 const int key_flag = cs_field_key_flag(key_id);
936 const int kr = cs_field_key_id_try("restart_name");
937
938 int n_diff = 0;
939
940 cs_lnum_t *old_key_val;
941 BFT_MALLOC(old_key_val, n_o_fields, cs_lnum_t);
942
943 char *sec_name;
944 BFT_MALLOC(sec_name, strlen("fields:") + strlen(key) + 1, char);
945 strcpy(sec_name, "fields:");
946 strcat(sec_name, key);
947
948 /* Read metadata */
949
950 retcode = cs_restart_check_section(r,
951 sec_name,
952 CS_MESH_LOCATION_NONE,
953 n_o_fields,
954 CS_TYPE_int);
955
956 if (retcode == CS_RESTART_SUCCESS)
957 retcode = cs_restart_read_section(r,
958 sec_name,
959 CS_MESH_LOCATION_NONE,
960 n_o_fields,
961 CS_TYPE_int,
962 old_key_val);
963
964 /* If data is available, compare models */
965
966 if (retcode == CS_RESTART_SUCCESS) {
967
968 for (int f_id = 0; f_id < n_fields; f_id++) {
969
970 const cs_field_t *f = cs_field_by_id(f_id);
971
972 if (key_flag == -1 || !(f->type & key_flag))
973 continue;
974
975 const char *f_name = NULL;
976 if (kr > -1)
977 f_name = cs_field_get_key_str(f, kr);
978 if (f_name == NULL)
979 f_name = f->name;
980
981 int old_f_id = cs_map_name_to_id_try(old_field_map, f_name);
982 if (old_f_id > -1) {
983 if (cs_field_get_key_int(f, key_id) != old_key_val[old_f_id])
984 n_diff += 1;
985 }
986
987 }
988
989 }
990 else if (retcode == CS_RESTART_ERR_EXISTS)
991 n_diff = -1;
992
993 else
994 bft_error
995 (__FILE__, __LINE__, 0,
996 _("Error %d reading \"%s\" in restart file \"%s\"."),
997 retcode, sec_name, cs_restart_get_name(r));
998
999 BFT_FREE(sec_name);
1000 BFT_FREE(old_key_val);
1001
1002 return n_diff;
1003 }
1004
1005 /*----------------------------------------------------------------------------
1006 * Check old turbulent flux model
1007 *
1008 * parameters:
1009 * r <-- associated restart file pointer
1010 * old_field_map <-- name to id map of fields in restart file
1011 *----------------------------------------------------------------------------*/
1012
1013 static void
_check_turb_flux_model(cs_restart_t * r,const cs_map_name_to_id_t * old_field_map)1014 _check_turb_flux_model(cs_restart_t *r,
1015 const cs_map_name_to_id_t *old_field_map)
1016 {
1017 const int n_fields = cs_field_n_fields();
1018
1019 int n_diff = _check_field_model(r,
1020 "turbulent_flux_model",
1021 old_field_map);
1022
1023 /* Read in legacy mode if required */
1024
1025 if (n_diff < 0) {
1026
1027 int kold = cs_field_key_id_try("old_scalar_num");
1028 const int key_id = cs_field_key_id("turbulent_flux_model");
1029
1030 if (kold > -1) {
1031
1032 n_diff = 0;
1033
1034 for (int f_id = 0; f_id < n_fields; f_id++) {
1035
1036 const cs_field_t *f = cs_field_by_id(f_id);
1037
1038 if (!(f->type & CS_FIELD_VARIABLE))
1039 continue;
1040
1041 int s_num = cs_field_get_key_int(f, kold);
1042
1043 if (s_num > 0) {
1044 char sec_name[128];
1045 cs_lnum_t old_s_model[1];
1046 snprintf(sec_name, 127, "turbulent_flux_model%04d", s_num);
1047 sec_name[127] = '\0';
1048 int retcode = cs_restart_read_section(r,
1049 sec_name,
1050 CS_MESH_LOCATION_NONE,
1051 1,
1052 CS_TYPE_int,
1053 old_s_model);
1054 if (retcode == CS_RESTART_SUCCESS) {
1055 if (cs_field_get_key_int(f, key_id) != old_s_model[0])
1056 n_diff += 1;
1057 }
1058 }
1059
1060 }
1061
1062 }
1063
1064 }
1065
1066 if (n_diff > 0)
1067 bft_printf
1068 (_("\n"
1069 " Warning: the turbulent flux model has been changed\n"
1070 " for %d fields relative to the restart file\n\n"
1071 " The computation continues, with a partial restart.\n"),
1072 n_diff);
1073 }
1074
1075 /*----------------------------------------------------------------------------
1076 * Read model option from file with compatibility for older files.
1077 *
1078 * parameters:
1079 * r <-- associated restart file pointer
1080 * m_name <-- name of model in current vesions of the restart file
1081 * m_name_old <-- name of model in different vesions of the restart file
1082 * count <-- number of values with current name
1083 * count_old <-- number of values with older name
1084 * options <-> model flags
1085 *
1086 * returns:
1087 * number of values read
1088 *----------------------------------------------------------------------------*/
1089
1090 static int
_read_model_option_compat(cs_restart_t * r,const char * m_name,const char * m_name_old,int count,int count_old,cs_lnum_t options[])1091 _read_model_option_compat(cs_restart_t *r,
1092 const char *m_name,
1093 const char *m_name_old,
1094 int count,
1095 int count_old,
1096 cs_lnum_t options[])
1097 {
1098 int retcount = 0;
1099
1100 int retcode = cs_restart_check_section(r,
1101 m_name,
1102 CS_MESH_LOCATION_NONE,
1103 1,
1104 CS_TYPE_int);
1105
1106 if (retcode == CS_RESTART_ERR_EXISTS) {
1107 retcode = cs_restart_check_section(r,
1108 m_name_old,
1109 CS_MESH_LOCATION_NONE,
1110 count_old,
1111 CS_TYPE_int);
1112 if (retcode == CS_RESTART_SUCCESS) {
1113 cs_restart_read_section(r,
1114 m_name_old,
1115 CS_MESH_LOCATION_NONE,
1116 count_old,
1117 CS_TYPE_int,
1118 options);
1119 if (retcode == CS_RESTART_SUCCESS)
1120 retcount = count_old;
1121 }
1122 }
1123 else {
1124 cs_restart_read_section(r,
1125 m_name,
1126 CS_MESH_LOCATION_NONE,
1127 count,
1128 CS_TYPE_int,
1129 options);
1130 if (retcode == CS_RESTART_SUCCESS)
1131 retcount = count;
1132 }
1133
1134 return retcount;
1135 }
1136
1137 /*----------------------------------------------------------------------------
1138 * Read turbulence variable section to a 1d array of cells.
1139 *
1140 * This wrapper function is intended mainly to allow calls with simplifed
1141 * arguments.
1142 *
1143 * parameters:
1144 * r <-- associated restart file pointer
1145 * base_name <-- base name of section to read
1146 * old_name <-- old name of section to read
1147 * t_id <-- associated time id
1148 * vals --> pointer to values
1149 *
1150 * returns:
1151 * 1 in case of error, 0 in case of success
1152 *----------------------------------------------------------------------------*/
1153
1154 static int
_read_turb_array_1d_compat(cs_restart_t * r,const char * base_name,const char * old_name,int t_id,cs_real_t * vals)1155 _read_turb_array_1d_compat(cs_restart_t *r,
1156 const char *base_name,
1157 const char *old_name,
1158 int t_id,
1159 cs_real_t *vals)
1160 {
1161 char sec_name[128];
1162 snprintf(sec_name, 127, "%s::vals::%01d", base_name, t_id);
1163 sec_name[127] = '\0';
1164 int retcode = CS_RESTART_SUCCESS;
1165 int retval = 0;
1166
1167 if (t_id == 0) {
1168 char old_sec_name[128];
1169 snprintf(old_sec_name, 127, "%s_ce_phase01", old_name);
1170 old_sec_name[127] = '\0';
1171 retcode = cs_restart_read_section_compat(r,
1172 sec_name,
1173 old_sec_name,
1174 CS_MESH_LOCATION_CELLS,
1175 1,
1176 CS_TYPE_cs_real_t,
1177 vals);
1178 }
1179 else
1180 retcode = cs_restart_read_section(r,
1181 sec_name,
1182 CS_MESH_LOCATION_CELLS,
1183 1,
1184 CS_TYPE_cs_real_t,
1185 vals);
1186 if (retcode != CS_RESTART_SUCCESS)
1187 retval = 1;
1188
1189 return retval;
1190 }
1191
1192 /*----------------------------------------------------------------------------
1193 * Compare current and previous turbulence model.
1194 *
1195 * parameters:
1196 * r <-- associated restart file pointer
1197 *
1198 * returns:
1199 * previous turbulence model
1200 *----------------------------------------------------------------------------*/
1201
1202 static int
_read_turbulence_model(cs_restart_t * r)1203 _read_turbulence_model(cs_restart_t *r)
1204 {
1205 cs_lnum_t iturb_old = 0;
1206
1207 /* Determine previous turbulence model */
1208
1209 int n_opt_vals = _read_model_option_compat(r,
1210 "turbulence_model",
1211 "modele_turbulence_phase01",
1212 1,
1213 1,
1214 &iturb_old);
1215
1216 if (n_opt_vals != 1)
1217 iturb_old = -999;
1218
1219 return iturb_old;
1220 }
1221
1222 /*----------------------------------------------------------------------------
1223 * Read turbulence variables converted from another model.
1224 *
1225 * This function should be called once we have read all "matching"
1226 * variables, as it assumes variables common to several turbulence models
1227 * (for example k between k-epsilon and k-omega, or epsilon between
1228 * k-epsilon and Rij-epsilon) have already been read.
1229 *
1230 * parameters:
1231 * r <-- associated restart file pointer
1232 * iturb_cur <-- current turbulence model number
1233 * iturb_old <-- previous turbulence model number
1234 * t_id <-- associated time id
1235 * read_flag <-> flag to track fields read, or NULL;
1236 * set to 1 for fields read, -1 for fields
1237 * failed to read (size: n_fields)
1238 *----------------------------------------------------------------------------*/
1239
1240 static void
_read_and_convert_turb_variables(cs_restart_t * r,int iturb_cur,int iturb_old,int t_id,int read_flag[])1241 _read_and_convert_turb_variables(cs_restart_t *r,
1242 int iturb_cur,
1243 int iturb_old,
1244 int t_id,
1245 int read_flag[])
1246 {
1247 int err_sum = 0, warn_sum = 0;
1248
1249 const int itytur_cur = iturb_cur/10;
1250 const int itytur_old = iturb_old/10;
1251
1252 /* If the turbulence model has not changed, nothing to do */
1253
1254 if (iturb_cur == iturb_old)
1255 return;
1256
1257 /* Warn user that turbulence model has changed */
1258
1259 bft_printf
1260 (_("\n"
1261 " Warning: the turbulence model has been changed\n"
1262 " relative to the restart file.\n\n"
1263 " current model: %d\n"
1264 " previous model: %d\n\n"
1265 " The computation continues, with a partial an/or adapted restart.\n"),
1266 iturb_cur, (int)iturb_old);
1267
1268 const cs_lnum_t n_cells = cs_glob_mesh->n_cells;
1269
1270 cs_real_t *v_tmp;
1271 BFT_MALLOC(v_tmp, n_cells, cs_real_t);
1272
1273 /* Now convert variables if needed. When we do not know how to
1274 deduce turbulence variables from available info, we ignore
1275 them, and keep the default initializations. */
1276
1277 int t_mask = (t_id == 0) ? 1 : 2 << (t_id-1);
1278
1279 if (itytur_cur == 2) { /* k-eps or v2f (phi-fbar or BL-v2/k) to k-eps */
1280
1281 /* if (itytur_old == 5),
1282 k and epsilon are common (so already read), and phi or f-bar
1283 is ignored, so nothing left to do here */
1284
1285 if (itytur_old == 3) { /* Rij to k-epsilon (epsilon already read) */
1286
1287 cs_real_t *v_k = CS_F_(k)->vals[t_id];
1288
1289 cs_real_6_t *rij;
1290 BFT_MALLOC(rij, n_cells, cs_real_6_t);
1291
1292 err_sum += _read_rij(r, CS_MESH_LOCATION_CELLS, 0, rij);
1293
1294 for (cs_lnum_t i = 0; i < n_cells; i++)
1295 v_k[i] = 0.5 * (rij[i][0] + rij[i][1] + rij[i][2]);
1296
1297 BFT_FREE(rij);
1298
1299 if (err_sum == 0)
1300 read_flag[CS_F_(k)->id] += t_mask;
1301
1302 }
1303 else if (iturb_old == 60) { /* k-omega to k-epsilon (k already read) */
1304
1305 cs_real_t *v_eps = CS_F_(eps)->vals[t_id];
1306 const cs_real_t *v_k = CS_F_(eps)->vals[t_id];
1307
1308 err_sum += _read_turb_array_1d_compat(r, "omega", "omega", t_id, v_eps);
1309
1310 /* Now transform omega to epsilon */
1311
1312 for (cs_lnum_t i = 0; i < n_cells; i++)
1313 v_eps[i] = cs_turb_cmu*v_eps[i]*v_k[i];
1314
1315 if (err_sum == 0)
1316 read_flag[CS_F_(eps)->id] += t_mask;
1317 }
1318
1319 }
1320 else if (itytur_cur == 3) { /* New computation is in Rij-epsilon */
1321
1322 if (itytur_old == 2 || iturb_old == 50) { /* k-epsilon or v2f
1323 (phi-fbar or BL-v2/k) to-> Rij
1324 (epsilon already read) */
1325
1326 cs_real_6_t *v_rij = (cs_real_6_t *)(CS_F_(rij)->vals[t_id]);
1327 cs_real_t *v_k;
1328 BFT_MALLOC(v_k, n_cells, cs_real_t);
1329
1330 err_sum += _read_turb_array_1d_compat(r, "k", "k", t_id, v_k);
1331
1332 double d2s3 = 2./3.;
1333
1334 for (cs_lnum_t i = 0; i < n_cells; i++) {
1335 double d2s3xk = v_k[i]*d2s3;
1336 v_rij[i][0] = d2s3xk;
1337 v_rij[i][1] = d2s3xk;
1338 v_rij[i][2] = d2s3xk;
1339 v_rij[i][3] = 0.;
1340 v_rij[i][4] = 0.;
1341 v_rij[i][5] = 0.;
1342 }
1343
1344 if (err_sum == 0)
1345 read_flag[CS_F_(rij)->id] += t_mask;
1346
1347 BFT_FREE(v_k);
1348
1349 }
1350
1351 /* if (itytur_old == 3) Rij to Rij variant; nothing to do
1352 * (when switching to EBRSM to another model,
1353 * keep alpha to default initialization) */
1354
1355 else if (iturb_old == 60) { /* k-omega to Rij */
1356
1357 cs_real_6_t *v_rij = (cs_real_6_t *)(CS_F_(rij)->vals[t_id]);
1358 cs_real_t *v_eps = CS_F_(eps)->vals[t_id];
1359
1360 cs_real_t *v_k;
1361 BFT_MALLOC(v_k, n_cells, cs_real_t);
1362
1363 err_sum += _read_turb_array_1d_compat(r, "k", "k", t_id, v_k);
1364 err_sum += _read_turb_array_1d_compat(r, "omega", "omega", t_id, v_eps);
1365
1366 /* Now transform omega to epsilon */
1367
1368 for (cs_lnum_t i = 0; i < n_cells; i++)
1369 v_eps[i] = cs_turb_cmu*v_eps[i]*v_k[i];
1370
1371 double d2s3 = 2./3.;
1372
1373 for (cs_lnum_t i = 0; i < n_cells; i++) {
1374 double d2s3xk = v_k[i]*d2s3;
1375 v_rij[i][0] = d2s3xk;
1376 v_rij[i][1] = d2s3xk;
1377 v_rij[i][2] = d2s3xk;
1378 v_rij[i][3] = 0.;
1379 v_rij[i][4] = 0.;
1380 v_rij[i][5] = 0.;
1381 }
1382
1383 if (err_sum == 0)
1384 read_flag[CS_F_(rij)->id] += t_mask;
1385
1386 BFT_FREE(v_k);
1387 }
1388
1389 } else if (itytur_cur == 5) { /* New computation is in v2f; */
1390
1391 /* if (itytur_old == 2) k-epsilon to v2f (phi-fbar or BL-v2/k)
1392 * k and epsilon already read, keep default
1393 * initializations for phi and f_bar or alpha,
1394 * so nothing special to do here */
1395
1396 if (itytur_old == 3) { /* Rij to v2f (phi-fbar or BL-v2/k)
1397 epsilon alread read; keep default initializations
1398 for phi and f_bar or alpha (v2 in phi is not a
1399 true Rij component) */
1400
1401 cs_real_t *v_k = CS_F_(k)->vals[t_id];
1402
1403 err_sum += _read_turb_array_1d_compat(r, "r11", "R11", t_id, v_k);
1404
1405 err_sum += _read_turb_array_1d_compat(r, "r22", "R22", t_id, v_tmp);
1406 for (cs_lnum_t i = 0; i < n_cells; i++)
1407 v_k[i] += v_tmp[i];
1408
1409 err_sum += _read_turb_array_1d_compat(r, "r33", "R33", t_id, v_tmp);
1410 for (cs_lnum_t i = 0; i < n_cells; i++)
1411 v_k[i] = 0.5 * (v_k[i] + v_tmp[i]);
1412
1413 if (err_sum == 0)
1414 read_flag[CS_F_(k)->id] += t_mask;
1415 }
1416
1417 /* if (itytur_old == 5) Switch between v2f variants; k, epsilon,
1418 * and phi already read; keep default initialization
1419 * for f_bar or alpha */
1420
1421 else if (iturb_old == 60) { /* Switch from k-omega to v2f;
1422 k already read; keep default initialization
1423 for f_bar or alpha */
1424
1425 cs_real_t *v_eps = CS_F_(eps)->vals[t_id];
1426 const cs_real_t *v_k = CS_F_(k)->vals[t_id];
1427
1428 err_sum += _read_turb_array_1d_compat(r, "omega", "omega", t_id, v_eps);
1429
1430 /* Now transform omega to epsilon */
1431
1432 for (cs_lnum_t i = 0; i < n_cells; i++)
1433 v_eps[i] = cs_turb_cmu*v_eps[i]*v_k[i];
1434
1435 if (err_sum == 0)
1436 read_flag[CS_F_(eps)->id] += t_mask;
1437 }
1438
1439 } else if (iturb_cur == 60) { /* New computation is in k-omega */
1440
1441 if (itytur_old == 2 || iturb_old == 50) { /* k-epsilon or v2f to k-omega;
1442 k already read */
1443
1444 cs_real_t *v_omg = CS_F_(omg)->vals[t_id];
1445 const cs_real_t *v_k = CS_F_(k)->vals[t_id];
1446
1447 err_sum += _read_turb_array_1d_compat(r, "epsilon", "eps", t_id, v_omg);
1448
1449 /* Now transform epsilon to omega */
1450
1451 for (cs_lnum_t i = 0; i < n_cells; i++)
1452 v_omg[i] /= (cs_turb_cmu*v_k[i]);
1453
1454 if (err_sum == 0)
1455 read_flag[CS_F_(omg)->id] += t_mask;
1456
1457 } else if (itytur_old == 3) { /* Rij to k-omega */
1458
1459 cs_real_t *v_k = CS_F_(k)->vals[t_id];
1460 cs_real_t *v_omg = CS_F_(omg)->vals[t_id];
1461
1462 cs_real_6_t *rij;
1463 BFT_MALLOC(rij, n_cells, cs_real_6_t);
1464
1465 err_sum += _read_rij(r, CS_MESH_LOCATION_CELLS, 0, rij);
1466
1467 for (cs_lnum_t i = 0; i < n_cells; i++)
1468 v_k[i] = 0.5 * (rij[i][0] + rij[i][1] + rij[i][2]);
1469
1470 BFT_FREE(rij);
1471
1472 err_sum += _read_turb_array_1d_compat(r, "epsilon", "eps", t_id, v_omg);
1473
1474 /* Now transform epsilon to omega */
1475
1476 for (cs_lnum_t i = 0; i < n_cells; i++)
1477 v_omg[i] /= (cs_turb_cmu*v_k[i]);
1478
1479 if (err_sum == 0) {
1480 read_flag[CS_F_(k)->id] += t_mask;
1481 read_flag[CS_F_(omg)->id] += t_mask;
1482 }
1483
1484 }
1485
1486 } else if (iturb_cur == 70) { /* current computation with the
1487 Spalart Allmaras (SA) model */
1488
1489 /* TODO perform the conversion from other models to SA. */
1490
1491 }
1492 else if (itytur_cur == 4) { /* LES mode */
1493
1494 if (itytur_old != 4) { /* restart from RANS */
1495
1496 const cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
1497 const cs_real_3_t *cell_cen = (const cs_real_3_t *)mq->cell_cen;
1498
1499 cs_real_3_t *v_vel = (cs_real_3_t *)(CS_F_(vel)->vals[t_id]);
1500
1501 cs_real_t *v_k;
1502 BFT_MALLOC(v_k, n_cells, cs_real_t);
1503 const cs_lnum_t n_cells_ext = cs_glob_mesh->n_cells_with_ghosts;
1504 bool use_previous_t = false;
1505 int inc = 1;
1506 cs_real_t *v_eps; /* Dissipation rate from the previous simulation */
1507 BFT_MALLOC(v_eps, n_cells, cs_real_t);
1508 cs_real_6_t *rst; /* Reynolds Stress Tensor for each cell */
1509 BFT_MALLOC(rst, n_cells , cs_real_6_t);
1510 cs_real_33_t *gradv; /* Velocity gradient for each cell */
1511 BFT_MALLOC(gradv, n_cells_ext, cs_real_33_t);
1512
1513 cs_field_gradient_vector(CS_F_(vel),
1514 use_previous_t,
1515 inc,
1516 gradv); /* Calculate velocity gradient */
1517
1518 warn_sum += _read_turb_array_1d_compat(r, "epsilon", "epsilon",
1519 t_id, v_eps);
1520 if (itytur_old == 3) { /* Rij */
1521
1522 warn_sum += _read_rij(r, CS_MESH_LOCATION_CELLS, t_id, rst);
1523
1524 }
1525
1526 /* Eddy viscosity model */
1527
1528 else {
1529
1530 warn_sum += _read_turb_array_1d_compat(r, "k", "k", t_id, v_k);
1531
1532 if (iturb_old == 60) { /* transform omega to epsilon */
1533 err_sum += _read_turb_array_1d_compat(r, "omega", "omega", t_id,
1534 v_eps);
1535
1536 for (cs_lnum_t i = 0; i < n_cells; i++)
1537 v_eps[i] = cs_turb_cmu*v_eps[i]*v_k[i];
1538 }
1539
1540 /* Loop over the cells to compute each component
1541 of the Reynolds stress tensor for each cell */
1542
1543 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
1544 cs_real_t divu = gradv[cell_id][0][0]
1545 + gradv[cell_id][1][1]
1546 + gradv[cell_id][2][2];
1547 // Turbulent viscosity = Experimental constant (0.09) * k^2 / epsilon
1548 cs_real_t nut = 0.09*cs_math_pow2(v_k[cell_id])/(v_eps[cell_id]);
1549
1550 // Diagonal of the Reynolds stress tensor
1551 cs_real_t xdiag = 2. / 3. *(v_k[cell_id]+ nut*divu);
1552
1553 rst[cell_id][0] = xdiag - 2.*nut*gradv[cell_id][0][0];
1554 rst[cell_id][1] = xdiag - 2.*nut*gradv[cell_id][1][1];
1555 rst[cell_id][2] = xdiag - 2.*nut*gradv[cell_id][2][2];
1556 rst[cell_id][3] = -nut*(gradv[cell_id][1][0]+gradv[cell_id][0][1]);
1557 rst[cell_id][4] = -nut*(gradv[cell_id][2][1]+gradv[cell_id][1][2]);
1558 rst[cell_id][5] = -nut*(gradv[cell_id][2][0]+gradv[cell_id][0][2]);
1559
1560 /* Clip it if necessary to get a SPD matrix
1561 * (same code as in clprij.f90) */
1562
1563 cs_real_t trrij = 2. * v_k[cell_id];
1564 /* Dimension less tensor R/2k */
1565 cs_real_t tensor[6];
1566 for (int i = 0; i < 6; i++)
1567 tensor[i] = rst[cell_id][i] / trrij;
1568
1569 cs_real_t eigen_vals[3];
1570 cs_math_sym_33_eigen(tensor, eigen_vals);
1571
1572 cs_real_t eigen_tol = 1.e-4;
1573 cs_real_t eigen_min = eigen_vals[0];
1574 cs_real_t eigen_max = eigen_vals[0];
1575 for (int i = 1; i < 3; i++) {
1576 eigen_min = CS_MIN(eigen_min, eigen_vals[i]);
1577 eigen_max = CS_MAX(eigen_max, eigen_vals[i]);
1578 }
1579
1580 /* If negative eigen value, return to isotropy */
1581 if ( eigen_min <= (eigen_tol*eigen_max)
1582 || eigen_min < cs_math_epzero) {
1583
1584 eigen_min = CS_MIN(eigen_min, - eigen_tol);
1585 cs_real_t eigen_offset
1586 = fmin(- eigen_min / (1./3. - eigen_min) + 0.1, 1.);
1587
1588 for (int i = 0; i < 6; i++) {
1589 rst[cell_id][i] *= (1. - eigen_offset);
1590
1591 /* Diagonal terms */
1592 if (i < 3)
1593 rst[cell_id][i] += trrij * (eigen_offset + eigen_tol) / 3.;
1594 }
1595 }
1596
1597 }
1598 }
1599
1600 /* Synthetic Eddy Method: eddies generation over the whole domain.
1601 Theory for the signal computation available in
1602 "A synthetic-eddy-method for generating inflow conditions
1603 for large-eddy simulations", Jarrin & al., 2006. */
1604
1605 /* FIXME: it is not clear whether everything is initialized;
1606 this code should maybe be moved to a function in cs_les_inflow.c,
1607 using the full inflow structure defined for the computation,
1608 or else use a specific function from that code with a
1609 simpler initialization path. */
1610
1611 int n_structures = cs_les_synthetic_eddy_get_n_restart_structures();
1612
1613 cs_inflow_sem_t *sem_in;
1614 BFT_MALLOC(sem_in, 1, cs_inflow_sem_t);
1615 sem_in->n_structures = n_structures;
1616 sem_in->volume_mode = 1;
1617 BFT_MALLOC(sem_in->position, sem_in->n_structures, cs_real_3_t);
1618 BFT_MALLOC(sem_in->energy, sem_in->n_structures, cs_real_3_t);
1619
1620 /* Velocity fluctuations before modifications with Lund's method */
1621 cs_real_3_t *fluctuations = NULL;
1622 BFT_MALLOC(fluctuations, n_cells, cs_real_3_t);
1623 cs_array_set_value_real(n_cells, 3, 0, (cs_real_t *)fluctuations);
1624
1625 cs_real_3_t *vel_l = NULL;
1626 BFT_MALLOC(vel_l, n_cells, cs_real_3_t);
1627 cs_array_set_value_real(n_cells, 3, 0, (cs_real_t *)vel_l);
1628
1629 cs_real_3_t *point_coordinates = NULL;
1630 BFT_MALLOC(point_coordinates, n_cells, cs_real_3_t);
1631 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
1632 for (cs_lnum_t j = 0; j < 3; j++)
1633 point_coordinates[cell_id][j] = cell_cen[cell_id][j];
1634 }
1635
1636 cs_real_t *point_weight = NULL;
1637 int initialize = 1;
1638 int verbosity = 1;
1639 cs_real_t t_cur = 0;
1640
1641 cs_les_synthetic_eddy_method(n_cells,
1642 NULL,
1643 point_coordinates,
1644 point_weight,
1645 initialize, verbosity,
1646 sem_in,
1647 t_cur,
1648 vel_l,
1649 rst,
1650 v_eps,
1651 fluctuations);
1652 cs_les_rescale_fluctuations(n_cells, rst, fluctuations);
1653
1654 /* Cancel fluctuations in solid zones */
1655 _cancel_in_solid_zones(3, (cs_real_t *)fluctuations);
1656
1657 for (cs_lnum_t cell_id = 0; cell_id < n_cells ; cell_id++) {
1658 /* Final update of velocities components unew = urans + u' */
1659 for (cs_lnum_t j = 0; j < 3; j++)
1660 v_vel[cell_id][j] += fluctuations[cell_id][j];
1661 }
1662
1663 BFT_FREE(fluctuations);
1664 BFT_FREE(vel_l);
1665 BFT_FREE(point_coordinates);
1666 BFT_FREE(rst);
1667 BFT_FREE(gradv);
1668 BFT_FREE(v_k);
1669 BFT_FREE(v_eps);
1670 BFT_FREE(sem_in->position);
1671 BFT_FREE(sem_in->energy);
1672 BFT_FREE(sem_in);
1673
1674 }
1675
1676 }
1677
1678 if (err_sum != 0)
1679 bft_error
1680 (__FILE__, __LINE__, 0,
1681 _("Error reading turbulence variables from previous model\n"
1682 "in restart file \"%s\"."),
1683 cs_restart_get_name(r));
1684
1685 BFT_FREE(v_tmp);
1686 }
1687
1688 /*============================================================================
1689 * Fortran wrapper function definitions
1690 *============================================================================*/
1691
1692 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
1693
1694 /*=============================================================================
1695 * Public function definitions
1696 *============================================================================*/
1697
1698 /*----------------------------------------------------------------------------*/
1699 /*!
1700 * \brief Read field metadata from checkpoint.
1701 *
1702 * Old ids associated to each field are determined for future use.
1703 * Note that when reading legacy files (Code_Saturne version 3.3 and below),
1704 * the old id will actually be the old scalar id (-1 for others).
1705 *
1706 * \param[in, out] r associated restart file pointer
1707 * \param[out] old_field_map name to id map of fields in restart file
1708 */
1709 /*----------------------------------------------------------------------------*/
1710
1711 void
cs_restart_read_field_info(cs_restart_t * r,cs_map_name_to_id_t ** old_field_map)1712 cs_restart_read_field_info(cs_restart_t *r,
1713 cs_map_name_to_id_t **old_field_map)
1714 {
1715 cs_lnum_t sizes[2];
1716 int retcode;
1717
1718 /* Initialization */
1719
1720 const int n_fields = cs_field_n_fields();
1721
1722 *old_field_map = NULL;
1723
1724 /* Now read field names, in id order */
1725
1726 int *type_buf = NULL;
1727 char *name_buf = NULL;
1728
1729 retcode = cs_restart_read_section(r,
1730 "fields:sizes",
1731 CS_MESH_LOCATION_NONE,
1732 2,
1733 CS_TYPE_int,
1734 sizes);
1735
1736 if (retcode == CS_RESTART_SUCCESS) {
1737
1738 const int n_o_fields = sizes[0];
1739 cs_lnum_t n_old[] = {0, 0, 0, 0};
1740 cs_lnum_t n_cur[] = {0, 0, 0, 0};
1741
1742 /* Now read main metadata */
1743
1744 BFT_MALLOC(name_buf, sizes[1] + 1, char);
1745 BFT_MALLOC(type_buf, sizes[0], int);
1746
1747 retcode = cs_restart_read_section(r,
1748 "fields:names",
1749 CS_MESH_LOCATION_NONE,
1750 sizes[1],
1751 CS_TYPE_char,
1752 name_buf);
1753
1754 if (retcode == CS_RESTART_SUCCESS) {
1755
1756 cs_map_name_to_id_t *o_map = cs_map_name_to_id_create();
1757
1758 const char *name = name_buf;
1759
1760 for (int i = 0, j = 0; j < n_o_fields; i++) {
1761 if (name_buf[i] == '\0') {
1762 cs_map_name_to_id(o_map, name);
1763 name = name_buf + i + 1;
1764 j++;
1765 }
1766 }
1767
1768 *old_field_map = o_map;
1769
1770 retcode = cs_restart_read_section(r,
1771 "fields:types",
1772 CS_MESH_LOCATION_NONE,
1773 sizes[0],
1774 CS_TYPE_int,
1775 type_buf);
1776
1777 if (retcode != CS_RESTART_SUCCESS) {
1778 for (int i = 0; i < n_fields; i++)
1779 type_buf[i] = 0;
1780 }
1781 }
1782
1783 BFT_FREE(name_buf);
1784
1785 /* Count variables */
1786
1787 for (int f_id = 0; f_id < n_fields; f_id++) {
1788 const cs_field_t *f = cs_field_by_id(f_id);
1789 if (f->type & CS_FIELD_VARIABLE) {
1790 if (f->type & CS_FIELD_USER)
1791 n_cur[1] += 1;
1792 else
1793 n_cur[0] += 1;
1794 }
1795 else if (f->type & CS_FIELD_PROPERTY) {
1796 if (f->type & CS_FIELD_USER)
1797 n_cur[3] += 1;
1798 else
1799 n_cur[2] += 1;
1800 }
1801 }
1802
1803 for (int f_id = 0; f_id < sizes[0]; f_id++) {
1804 if (type_buf[f_id] & CS_FIELD_VARIABLE) {
1805 if (type_buf[f_id] & CS_FIELD_USER)
1806 n_old[1] += 1;
1807 else
1808 n_old[0] += 1;
1809 }
1810 else if (type_buf[f_id] & CS_FIELD_PROPERTY) {
1811 if (type_buf[f_id] & CS_FIELD_USER)
1812 n_old[3] += 1;
1813 else
1814 n_old[2] += 1;
1815 }
1816 }
1817
1818 /* Warn in case of change */
1819
1820 if ( n_cur[0] != n_old[0] || n_cur[1] != n_old[1]
1821 || n_cur[2] != n_old[2] || n_cur[3] != n_old[3])
1822 bft_printf
1823 (_("\n"
1824 " Warning: the number of variables or properties has been changed\n"
1825 " relative to the restart file.\n\n"
1826 " current: %d model and %d user variables\n"
1827 " %d model and %d user properties\n"
1828 " previous: %d model and %d user variables\n"
1829 " %d model and %d user properties\n\n"
1830 " The computation continues, with a partial restart.\n"),
1831 (int)(n_cur[0]), (int)(n_cur[1]), (int)(n_cur[2]), (int)(n_cur[3]),
1832 (int)(n_old[0]), (int)(n_old[1]), (int)(n_old[2]), (int)(n_old[3]));
1833
1834 BFT_FREE(type_buf);
1835 }
1836
1837 /* Read legacy metadata */
1838
1839 else if (cs_restart_is_from_ncfd() == 0)
1840 _read_legacy_field_info(r);
1841 }
1842
1843 /*----------------------------------------------------------------------------*/
1844 /*!
1845 * \brief Write field metadata to main checkpoint.
1846 *
1847 * \param[in, out] r associated restart file pointer
1848 */
1849 /*----------------------------------------------------------------------------*/
1850
1851 void
cs_restart_write_field_info(cs_restart_t * r)1852 cs_restart_write_field_info(cs_restart_t *r)
1853 {
1854 cs_lnum_t sizes[2];
1855
1856 int n_fields = cs_field_n_fields();
1857
1858 sizes[0] = n_fields;
1859 sizes[1] = 0;
1860
1861 /* Build field name buffer */
1862
1863 for (int f_id = 0; f_id < n_fields; f_id++) {
1864 const cs_field_t *f = cs_field_by_id(f_id);
1865 sizes[1] += strlen(f->name) + 1;
1866 }
1867
1868 cs_lnum_t *type_buf;
1869 char *name_buf;
1870
1871 BFT_MALLOC(type_buf, n_fields, cs_lnum_t);
1872 BFT_MALLOC(name_buf, sizes[1] + 1, char);
1873
1874 sizes[1] = 0;
1875
1876 for (int f_id = 0; f_id < n_fields; f_id++) {
1877
1878 const cs_field_t *f = cs_field_by_id(f_id);
1879
1880 size_t l = strlen(f->name) + 1;
1881 memcpy(name_buf + sizes[1], f->name, l);
1882 sizes[1] += l;
1883
1884 type_buf[f_id] = f->type;
1885
1886 }
1887
1888 /* Now write data */
1889
1890 cs_restart_write_section(r,
1891 "fields:sizes",
1892 CS_MESH_LOCATION_NONE,
1893 2,
1894 CS_TYPE_int,
1895 sizes);
1896
1897 cs_restart_write_section(r,
1898 "fields:names",
1899 CS_MESH_LOCATION_NONE,
1900 sizes[1],
1901 CS_TYPE_char,
1902 name_buf);
1903
1904 cs_restart_write_section(r,
1905 "fields:types",
1906 CS_MESH_LOCATION_NONE,
1907 n_fields,
1908 CS_TYPE_int,
1909 type_buf);
1910
1911 BFT_FREE(name_buf);
1912 BFT_FREE(type_buf);
1913
1914 bft_printf(_(" Wrote field names and types to checkpoint"
1915 " at iteration %d: %s\n"),
1916 cs_glob_time_step->nt_cur, cs_restart_get_name(r));
1917 }
1918
1919 /*----------------------------------------------------------------------------*/
1920 /*!
1921 * \brief Read variables from checkpoint.
1922 *
1923 * \param[in, out] r associated restart file pointer
1924 * \param[in] old_field_map name to id map of fields in restart file
1925 * \param[in] t_id_flag -1: all time values; 0: current values;
1926 * > 0: previous values
1927 * \param[in, out] read_flag optional flag to track fields read,
1928 * or NULL; set to sum of 2^time_id for fields
1929 * read (size: n_fields)
1930 */
1931 /*----------------------------------------------------------------------------*/
1932
1933 void
cs_restart_read_variables(cs_restart_t * r,const cs_map_name_to_id_t * old_field_map,int t_id_flag,int read_flag[])1934 cs_restart_read_variables(cs_restart_t *r,
1935 const cs_map_name_to_id_t *old_field_map,
1936 int t_id_flag,
1937 int read_flag[])
1938 {
1939 const int n_fields = cs_field_n_fields();
1940
1941 /* Initialization */
1942
1943 int *_read_flag = read_flag;
1944
1945 if (_read_flag == NULL) {
1946 BFT_MALLOC(_read_flag, n_fields, int);
1947 for (int f_id = 0; f_id < n_fields; f_id++)
1948 _read_flag[f_id] = 0;
1949 }
1950
1951 /* Read metadata for turbulent flux models
1952 (turbulent flux model useful only for warnings, as the turbulent flux
1953 field id is sufficient to determine data to be read) */
1954
1955 _check_turb_flux_model(r, old_field_map);
1956
1957 /* Read field data */
1958
1959 for (int f_id = 0; f_id < n_fields; f_id++) {
1960
1961 /* Base fields */
1962
1963 const cs_field_t *f = cs_field_by_id(f_id);
1964 if (f->type & CS_FIELD_VARIABLE) {
1965 int t_id_s = 0;
1966 int t_id_e = f->n_time_vals;
1967 if (t_id_flag == 0)
1968 t_id_e = 1;
1969 else if (t_id_flag > 0)
1970 t_id_s = 1;
1971
1972 for (int t_id = t_id_s; t_id < t_id_e; t_id++) {
1973 int t_mask = (t_id == 0) ? 1 : 2 << (t_id-1);
1974 if (_read_flag[f->id] & t_mask)
1975 continue;
1976
1977 int retval = cs_restart_read_field_vals(r, f_id, t_id);
1978
1979 if (retval == CS_RESTART_SUCCESS)
1980 _read_flag[f_id] += t_mask;
1981 }
1982
1983 }
1984 }
1985
1986 /* Read linked field data;
1987 note that when turbulent fluxes are variables and not properties,
1988 they should already have been read by the main loop on variables,
1989 so use of _read_flag here should avoid writing twice.
1990 */
1991
1992 cs_restart_read_linked_fields(r,
1993 old_field_map,
1994 "turbulent_flux_id",
1995 _read_flag);
1996
1997 cs_restart_read_linked_fields(r,
1998 old_field_map,
1999 "gwf_sorbed_concentration_id",
2000 _read_flag);
2001
2002 cs_restart_read_linked_fields(r,
2003 old_field_map,
2004 "gwf_precip_concentration_id",
2005 _read_flag);
2006
2007 /* Read and convert turbulence variables in case of model change */
2008
2009 const cs_turb_model_t *turb_model = cs_get_glob_turb_model();
2010 assert(turb_model != NULL);
2011 const int iturb_cur = turb_model->iturb;
2012 const int iturb_old = _read_turbulence_model(r);
2013
2014 if (iturb_cur != iturb_old)
2015 _read_and_convert_turb_variables(r, iturb_cur, iturb_old, 0, _read_flag);
2016
2017 /* Warning or error for main unread variables */
2018
2019 if (t_id_flag < 1) {
2020
2021 if (cs_glob_field_pointers != NULL) {
2022
2023 /* standard calculation */
2024 if (CS_F_(p)) {
2025 if (!(_read_flag[CS_F_(vel)->id] & 1) || !(_read_flag[CS_F_(p)->id] & 1))
2026 bft_error
2027 (__FILE__, __LINE__, 0,
2028 _("Error reading velocity/pressure values in restart file \"%s\"."),
2029 cs_restart_get_name(r));
2030 /* ground water flow calculation */
2031 } else if (CS_F_(head)) {
2032 if ( !(_read_flag[CS_F_(vel)->id] & 1)
2033 || !(_read_flag[CS_F_(head)->id] & 1))
2034 bft_error
2035 (__FILE__, __LINE__, 0,
2036 _("Error reading velocity/hydraulic head values in restart "
2037 "file \"%s\"."),
2038 cs_restart_get_name(r));
2039 }
2040
2041 }
2042
2043 int n_missing = 0;
2044
2045 for (int f_id = 0; f_id < n_fields; f_id++) {
2046 const cs_field_t *f = cs_field_by_id(f_id);
2047 if (f->type & CS_FIELD_VARIABLE) {
2048 if (!(_read_flag[f_id] & 1))
2049 n_missing++;
2050 }
2051 }
2052
2053 if (n_missing > 0) {
2054 bft_printf
2055 (_("\n"
2056 " Warning: the following variables could not be found or read\n"
2057 " in restart file \"%s\", so default initializations\n"
2058 " will be used:\n\n"), cs_restart_get_name(r));
2059 for (int f_id = 0; f_id < n_fields; f_id++) {
2060 const cs_field_t *f = cs_field_by_id(f_id);
2061 if (f->type & CS_FIELD_VARIABLE) {
2062 if (!(_read_flag[f_id] & 1))
2063 bft_printf(" %s\n", cs_field_get_label(f));
2064 }
2065 }
2066 bft_printf("\n");
2067 }
2068
2069 }
2070
2071 /* Cleanup */
2072
2073 if (_read_flag != read_flag)
2074 BFT_FREE(_read_flag);
2075
2076 bft_printf(_(" Read variables from restart: %s\n"),
2077 cs_restart_get_name(r));
2078 }
2079
2080 /*----------------------------------------------------------------------------*/
2081 /*!
2082 * \brief Write variables to checkpoint.
2083 *
2084 * \param[in, out] r associated restart file pointer
2085 * \param[in] t_id_flag -1: all time values; 0: current values;
2086 * > 0: previous values
2087 * \param[in, out] write_flag optional flag to track fields written, or NULL;
2088 * set to sum of 2^time_id for fields written
2089 * (size: n_fields)
2090 */
2091 /*----------------------------------------------------------------------------*/
2092
2093 void
cs_restart_write_variables(cs_restart_t * r,int t_id_flag,int write_flag[])2094 cs_restart_write_variables(cs_restart_t *r,
2095 int t_id_flag,
2096 int write_flag[])
2097 {
2098 int n_fields = cs_field_n_fields();
2099
2100 /* Initialization */
2101
2102 int *_write_flag = write_flag;
2103
2104 if (_write_flag == NULL) {
2105 BFT_MALLOC(_write_flag, n_fields, int);
2106 for (int f_id = 0; f_id < n_fields; f_id++)
2107 _write_flag[f_id] = 0;
2108 }
2109
2110 /* Write metadata for turbulent flux models
2111 (turbulent flux model useful only for warnings, as the turbulent flux
2112 field id is sufficient to determine data to be written;
2113 last test on n_turbt is a minor optimization to possibly avoid a call
2114 to cs_restart_write_linked_fields, but that call would be no
2115 more costly than the loop to determine the turbulent flux model
2116 in the first case...) */
2117
2118 int n_turbt = 0;
2119
2120 {
2121 cs_lnum_t *turbt_buf;
2122
2123 BFT_MALLOC(turbt_buf, n_fields, cs_lnum_t);
2124
2125 for (int f_id = 0; f_id < n_fields; f_id++)
2126 turbt_buf[f_id] = 0;
2127
2128 int k_sca = cs_field_key_id("scalar_id");
2129 int k_turbt = cs_field_key_id("turbulent_flux_model");
2130
2131 for (int f_id = 0; f_id < n_fields; f_id++) {
2132 const cs_field_t *f = cs_field_by_id(f_id);
2133 if (f->type & CS_FIELD_VARIABLE) {
2134 int s_num = cs_field_get_key_int(f, k_sca);
2135 if (s_num > 0) {
2136 int f_turbt = cs_field_get_key_int(f, k_turbt);
2137 if (f_turbt > 0) {
2138 turbt_buf[f_id] = f_turbt;
2139 n_turbt++;
2140 }
2141 }
2142 }
2143 }
2144
2145 if (n_turbt > 0 && t_id_flag < 1)
2146 cs_restart_write_section(r,
2147 "fields:turbulent_flux_model",
2148 CS_MESH_LOCATION_NONE,
2149 n_fields,
2150 CS_TYPE_int,
2151 turbt_buf);
2152
2153 BFT_FREE(turbt_buf);
2154 }
2155
2156 /* Write field data */
2157
2158 for (int f_id = 0; f_id < n_fields; f_id++) {
2159
2160 /* Base fields */
2161
2162 const cs_field_t *f = cs_field_by_id(f_id);
2163 if (f->type & CS_FIELD_VARIABLE) {
2164 int t_id_s = 0;
2165 int t_id_e = f->n_time_vals;
2166 if (t_id_flag == 0)
2167 t_id_e = 1;
2168 else if (t_id_flag > 0)
2169 t_id_s = 1;
2170
2171 for (int t_id = t_id_s; t_id < t_id_e; t_id++) {
2172 int t_mask = (t_id == 0) ? 1 : 2 << (t_id-1);
2173 if (_write_flag[f_id] & t_mask)
2174 continue;
2175
2176 cs_restart_write_field_vals(r, f_id, t_id);
2177
2178 _write_flag[f_id] += t_mask;
2179 }
2180
2181 }
2182 }
2183
2184 /* Write linked field data;
2185 note that when turbulent fluxes are variables and not properties,
2186 they should already have been written by the main loop on variables,
2187 so use of _write_flag here should avoid writing twice.
2188 */
2189
2190 if (n_turbt > 0)
2191 cs_restart_write_linked_fields(r,
2192 "turbulent_flux_id",
2193 _write_flag);
2194
2195 cs_restart_write_linked_fields(r,
2196 "gwf_sorbed_concentration_id",
2197 _write_flag);
2198
2199 cs_restart_write_linked_fields(r,
2200 "gwf_precip_concentration_id",
2201 _write_flag);
2202
2203 /* Cleanup */
2204
2205 if (_write_flag != write_flag)
2206 BFT_FREE(_write_flag);
2207
2208 bft_printf(_(" Wrote main variables to checkpoint: %s\n"),
2209 cs_restart_get_name(r));
2210 }
2211
2212 /*----------------------------------------------------------------------------*/
2213 /*!
2214 * \brief Read fields depending on others from checkpoint.
2215 *
2216 * Old ids associate to each field are determined for future use.
2217 * Note that when reading legacy files (Code_Saturne version 3.3 and below),
2218 * the old id will actually be the old scalar id (-1 for others).
2219 *
2220 * \param[in, out] r associated restart file pointer
2221 * \param[in] old_field_map name to id map of fields in restart file
2222 * \param[in] key key for field association
2223 * \param[in, out] read_flag optional flag to track fields read, or NULL;
2224 * set to sum of 2^time_id for fields read,
2225 * -1 for fields failed to read (size: n_fields)
2226 */
2227 /*----------------------------------------------------------------------------*/
2228
2229 void
cs_restart_read_linked_fields(cs_restart_t * r,const cs_map_name_to_id_t * old_field_map,const char * key,int read_flag[])2230 cs_restart_read_linked_fields(cs_restart_t *r,
2231 const cs_map_name_to_id_t *old_field_map,
2232 const char *key,
2233 int read_flag[])
2234 {
2235 int retcode;
2236
2237 /* Initialization */
2238
2239 int n_required = 0, n_legacy_read = 0;
2240
2241 const int n_fields = cs_field_n_fields();
2242 const int n_o_fields = cs_map_name_to_id_size(old_field_map);
2243
2244 const int key_id = cs_field_key_id_try(key);
2245 const int key_flag = cs_field_key_flag(key_id);
2246 const int kr = cs_field_key_id_try("restart_name");
2247
2248 /* First, check if we need to read anything */
2249
2250 for (int f_id = 0; f_id < n_fields; f_id++) {
2251 const cs_field_t *f = cs_field_by_id(f_id);
2252 if (key_flag != 0) {
2253 if (key_flag == -1 || !(f->type & key_flag))
2254 continue;
2255 }
2256 if (cs_field_get_key_int(f, key_id) > -1)
2257 n_required += 1;
2258 }
2259
2260 if (n_required < 1)
2261 return;
2262
2263 /* Now try to read */
2264
2265 int *_read_flag = read_flag;
2266
2267 if (_read_flag == NULL) {
2268 BFT_MALLOC(_read_flag, n_fields, int);
2269 for (int f_id = 0; f_id < n_fields; f_id++)
2270 _read_flag[f_id] = 0;
2271 }
2272
2273 cs_lnum_t *old_key_val;
2274 BFT_MALLOC(old_key_val, n_o_fields, cs_lnum_t);
2275
2276 char *sec_name;
2277 BFT_MALLOC(sec_name, strlen("fields:") + strlen(key) + 1, char);
2278 strcpy(sec_name, "fields:");
2279 strcat(sec_name, key);
2280
2281 /* Read metadata */
2282
2283 retcode = cs_restart_check_section(r,
2284 sec_name,
2285 CS_MESH_LOCATION_NONE,
2286 n_o_fields,
2287 CS_TYPE_int);
2288
2289 /* Try to read in compatibility mode if section not found */
2290
2291 if (retcode == CS_RESTART_ERR_EXISTS)
2292 n_legacy_read = _read_linked_fields_legacy(r, key, _read_flag);
2293
2294 /* Otherwise try to read section anyways to log warning */
2295
2296 if (n_legacy_read == 0)
2297 retcode = cs_restart_read_section(r,
2298 sec_name,
2299 CS_MESH_LOCATION_NONE,
2300 n_o_fields,
2301 CS_TYPE_int,
2302 old_key_val);
2303
2304 BFT_FREE(sec_name);
2305
2306 if (retcode == CS_RESTART_SUCCESS && n_legacy_read == 0) {
2307
2308 for (int f_id = 0; f_id < n_fields; f_id++) {
2309
2310 const cs_field_t *f = cs_field_by_id(f_id);
2311
2312 if (key_flag != 0) {
2313 if (key_flag == -1 || !(f->type & key_flag))
2314 continue;
2315 }
2316
2317 const int lnk_f_id = cs_field_get_key_int(f, key_id);
2318
2319 if (lnk_f_id > -1) {
2320
2321 cs_field_t *f_lnk = cs_field_by_id(lnk_f_id);
2322 const char *f_lnk_name = NULL;
2323
2324 if (_read_flag[lnk_f_id] != 0)
2325 continue;
2326
2327 /* Check if dependent field has explicit rename */
2328 if (kr > -1)
2329 f_lnk_name = cs_field_get_key_str(f, kr);
2330
2331 /* Otherwise, determine matching name,
2332 with possible parent field rename */
2333
2334 if (f_lnk_name == NULL) {
2335 const char *f_name = NULL;
2336 if (kr > -1)
2337 f_name = cs_field_get_key_str(f, kr);
2338 if (f_name == NULL)
2339 f_name = f->name;
2340 int old_f_id = cs_map_name_to_id_try(old_field_map, f_name);
2341 if (old_f_id > -1) {
2342 int old_lnk_f_id = old_key_val[old_f_id];
2343 if (old_lnk_f_id > -1)
2344 f_lnk_name
2345 = cs_map_name_to_id_reverse(old_field_map, old_lnk_f_id);
2346 else
2347 f_lnk_name = f_lnk->name;
2348 }
2349 }
2350
2351 /* Now we know which field name to read */
2352
2353 if (f_lnk_name != NULL) {
2354
2355 _read_flag[lnk_f_id] = -1;
2356
2357 for (int t_id = 0; t_id < f_lnk->n_time_vals; t_id++) {
2358 retcode = _read_field_vals(r, f_lnk_name, t_id, f_lnk);
2359 if (retcode == CS_RESTART_SUCCESS) {
2360 if (t_id == 0)
2361 _read_flag[lnk_f_id] = 1;
2362 else
2363 _read_flag[lnk_f_id] += (2 << (t_id-1));
2364 }
2365 else
2366 break;
2367 }
2368
2369 }
2370 else if (_read_flag[lnk_f_id] == 0) {
2371 _read_flag[lnk_f_id] = -1;
2372 bft_printf(_(" %s: no matching data for field \"%s\"\n"),
2373 cs_restart_get_name(r), f_lnk->name);
2374 if (f_lnk_name != NULL && f_lnk_name != f_lnk->name)
2375 bft_printf(_(" (was named \"%s\")\n"),
2376 f_lnk_name);
2377 }
2378
2379 }
2380
2381 }
2382 }
2383
2384 BFT_FREE(old_key_val);
2385
2386 if (read_flag != _read_flag)
2387 BFT_FREE(_read_flag);
2388 }
2389
2390 /*----------------------------------------------------------------------------*/
2391 /*!
2392 * \brief Write fields depending on others to checkpoint.
2393 *
2394 * \param[in, out] r associated restart file pointer
2395 * \param[in] key key for field association
2396 * \param[in, out] write_flag optional flag to track fields written, or NULL;
2397 * set to sum of 2^time_id for fields written
2398 * (size: n_fields)
2399 * \return number of fields written
2400 */
2401 /*----------------------------------------------------------------------------*/
2402
2403 int
cs_restart_write_linked_fields(cs_restart_t * r,const char * key,int write_flag[])2404 cs_restart_write_linked_fields(cs_restart_t *r,
2405 const char *key,
2406 int write_flag[])
2407 {
2408 int retval = 0;
2409
2410 /* Initialization */
2411
2412 const int n_fields = cs_field_n_fields();
2413
2414 const int key_id = cs_field_key_id_try(key);
2415 const int key_flag = cs_field_key_flag(key_id);
2416
2417 int *_write_flag = write_flag;
2418
2419 if (_write_flag == NULL) {
2420 BFT_MALLOC(_write_flag, n_fields, int);
2421 for (int f_id = 0; f_id < n_fields; f_id++)
2422 _write_flag[f_id] = 0;
2423 }
2424
2425 cs_lnum_t *key_val;
2426 BFT_MALLOC(key_val, n_fields, cs_lnum_t);
2427
2428 char *sec_name;
2429 BFT_MALLOC(sec_name, strlen("fields:") + strlen(key) + 1, char);
2430 strcpy(sec_name, "fields:");
2431 strcat(sec_name, key);
2432
2433 /* Write metadata */
2434
2435 for (int f_id = 0; f_id < n_fields; f_id++) {
2436 key_val[f_id] = -1;
2437 const cs_field_t *f = cs_field_by_id(f_id);
2438 if (key_flag != 0) {
2439 if (key_flag == -1 || !(f->type & key_flag))
2440 continue;
2441 }
2442 key_val[f_id] = cs_field_get_key_int(f, key_id);
2443 }
2444
2445 cs_restart_write_section(r,
2446 sec_name,
2447 CS_MESH_LOCATION_NONE,
2448 n_fields,
2449 CS_TYPE_int,
2450 key_val);
2451
2452 BFT_FREE(sec_name);
2453
2454 for (int f_id = 0; f_id < n_fields; f_id++) {
2455
2456 const int lnk_f_id = key_val[f_id];
2457 if (lnk_f_id < 0 || _write_flag[lnk_f_id] != 0)
2458 continue;
2459
2460 cs_field_t *f_lnk = cs_field_by_id(lnk_f_id);
2461
2462 _write_flag[lnk_f_id] = -1;
2463
2464 /* Now write field values */
2465
2466 for (int t_id = 0; t_id < f_lnk->n_time_vals; t_id++) {
2467
2468 cs_restart_write_field_vals(r, lnk_f_id, t_id);
2469
2470 if (t_id == 0)
2471 _write_flag[lnk_f_id] = 1;
2472 else
2473 _write_flag[lnk_f_id] += (2 << (t_id-1));
2474
2475 }
2476
2477 retval += 1;
2478
2479 }
2480
2481 BFT_FREE(key_val);
2482
2483 if (_write_flag != write_flag)
2484 BFT_FREE(_write_flag);
2485
2486 return retval;
2487 }
2488
2489 /*----------------------------------------------------------------------------*/
2490 /*!
2491 * \brief Read boundary condition coefficients for all fields from checkpoint.
2492 *
2493 * \param[in, out] r associated restart file pointer
2494 */
2495 /*----------------------------------------------------------------------------*/
2496
2497 void
cs_restart_read_bc_coeffs(cs_restart_t * r)2498 cs_restart_read_bc_coeffs(cs_restart_t *r)
2499 {
2500 int c_id, f_id;
2501
2502 int errcount = 0;
2503 const int coupled_key_id = cs_field_key_id_try("coupled");
2504 const int n_fields = cs_field_n_fields();
2505 char old_name_xx[128] = "", old_name_yy[128] = "", old_name_zz[128] = "";
2506 char old_name_xy[128] = "", old_name_yz[128] = "", old_name_xz[128] = "";
2507 const int kr = cs_field_key_id_try("restart_name");
2508
2509 /* Loop on all fields, to search for those defined on all cells
2510 and with BC coefficients */
2511
2512 for (f_id = 0; f_id < n_fields; f_id++) {
2513
2514 const cs_field_t *f = cs_field_by_id(f_id);
2515
2516 if ( f->location_id == CS_MESH_LOCATION_CELLS
2517 && f->bc_coeffs != NULL) {
2518
2519 /* Check for presence of coefficients */
2520
2521 int coupled = 0;
2522 int n_loc_vals = 1;
2523
2524 int32_t coeff_p[] = {0, 0, 0, 0, 0, 0, 0, 0};
2525
2526 cs_real_t *p[] = {f->bc_coeffs->a,
2527 f->bc_coeffs->b,
2528 f->bc_coeffs->af,
2529 f->bc_coeffs->bf,
2530 f->bc_coeffs->ad,
2531 f->bc_coeffs->bd,
2532 f->bc_coeffs->ac,
2533 f->bc_coeffs->bc};
2534
2535 for (c_id = 0; c_id < 8; c_id++) {
2536 if (p[c_id] != NULL) {
2537 coeff_p[c_id] = 1;
2538 /* avoid double reads/writes in case of aliasing */
2539 for (int i = 0; i < c_id; i++) {
2540 if (p[i] == p[c_id])
2541 coeff_p[c_id] = 0;
2542 }
2543 }
2544 }
2545
2546 cs_parall_max(8, CS_INT32, coeff_p);
2547
2548 if (f->dim > 1 && coupled_key_id > -1)
2549 coupled = cs_field_get_key_int(f, coupled_key_id);
2550
2551 for (c_id = 0; c_id < 8; c_id++) {
2552 int retval;
2553 char *sec_name = NULL;
2554 cs_real_t *c = p[c_id];
2555 const char *name = NULL;
2556 if (kr > -1)
2557 name = cs_field_get_key_str(f, kr);
2558 if (name == NULL)
2559 name = f->name;
2560 if (coeff_p[c_id] == 0)
2561 continue;
2562
2563 if (coupled) {
2564 if (c_id %2 == 0)
2565 n_loc_vals = f->dim;
2566 else
2567 n_loc_vals = f->dim * f->dim;
2568 }
2569 else { /* uncoupled */
2570 n_loc_vals = f->dim;
2571 }
2572
2573 BFT_MALLOC(sec_name,
2574 strlen(name) + strlen(_coeff_name[c_id]) + 3,
2575 char);
2576 sprintf(sec_name, "%s::%s", name, _coeff_name[c_id]);
2577
2578 retval = cs_restart_check_section(r,
2579 sec_name,
2580 f->location_id,
2581 f->dim,
2582 CS_TYPE_cs_real_t);
2583
2584 if (f->dim == 6 && retval == CS_RESTART_ERR_EXISTS) {
2585 sprintf(sec_name, "rij::%s", _coeff_name[c_id]);
2586 snprintf(old_name_xx, 127, "r11::%s", _coeff_name[c_id]);
2587 snprintf(old_name_yy, 127, "r22::%s", _coeff_name[c_id]);
2588 snprintf(old_name_zz, 127, "r33::%s", _coeff_name[c_id]);
2589 snprintf(old_name_xy, 127, "r12::%s", _coeff_name[c_id]);
2590 snprintf(old_name_yz, 127, "r23::%s", _coeff_name[c_id]);
2591 snprintf(old_name_xz, 127, "r13::%s", _coeff_name[c_id]);
2592 if (c_id %2 == 0) {
2593 retval = cs_restart_read_real_6_t_compat(r,
2594 sec_name,
2595 old_name_xx,
2596 old_name_yy,
2597 old_name_zz,
2598 old_name_xy,
2599 old_name_yz,
2600 old_name_xz,
2601 f->location_id,
2602 (cs_real_6_t *)(f->val));
2603 }
2604 else {
2605 retval = cs_restart_read_real_66_t_compat(r,
2606 sec_name,
2607 old_name_xx,
2608 old_name_yy,
2609 old_name_zz,
2610 old_name_xy,
2611 old_name_yz,
2612 old_name_xz,
2613 f->location_id,
2614 (cs_real_66_t *)(f->val));
2615 }
2616 }
2617 else {
2618 retval = cs_restart_read_section(r,
2619 sec_name,
2620 3, /* location_id */
2621 n_loc_vals,
2622 CS_TYPE_cs_real_t,
2623 c);
2624 }
2625 if (retval != CS_RESTART_SUCCESS)
2626 errcount += 1;
2627
2628 BFT_FREE(sec_name);
2629
2630 } /* End of loop in i (coeff type) */
2631
2632 } /* End for field with BC coeffs */
2633
2634 } /* End of loop on fields */
2635
2636 if (errcount > 0) {
2637 cs_base_warn(__FILE__, __LINE__);
2638 bft_printf(_("\nSome boundary condition coefficients "
2639 "could not be read from a restart file;\n"
2640 "they will be initialized with default values.\n\n"));
2641 }
2642
2643 }
2644
2645 /*----------------------------------------------------------------------------*/
2646 /*!
2647 * \brief Write boundary condition coefficients for all fields to checkpoint.
2648 *
2649 * \param[in, out] r associated restart file pointer
2650 */
2651 /*----------------------------------------------------------------------------*/
2652
2653 void
cs_restart_write_bc_coeffs(cs_restart_t * r)2654 cs_restart_write_bc_coeffs(cs_restart_t *r)
2655 {
2656 int c_id, f_id;
2657
2658 const int coupled_key_id = cs_field_key_id_try("coupled");
2659 const int n_fields = cs_field_n_fields();
2660
2661 /* Loop on all fields, to search for those defined on all cells
2662 and with BC coefficients */
2663
2664 for (f_id = 0; f_id < n_fields; f_id++) {
2665
2666 const cs_field_t *f = cs_field_by_id(f_id);
2667
2668 if ( f->location_id == CS_MESH_LOCATION_CELLS
2669 && f->bc_coeffs != NULL) {
2670
2671 /* Check for presence of coefficients */
2672
2673 int coupled = 0;
2674 int n_loc_vals = 1;
2675
2676 int32_t coeff_p[] = {0, 0, 0, 0, 0, 0, 0, 0};
2677 cs_real_t *p[] = {f->bc_coeffs->a,
2678 f->bc_coeffs->b,
2679 f->bc_coeffs->af,
2680 f->bc_coeffs->bf,
2681 f->bc_coeffs->ad,
2682 f->bc_coeffs->bd,
2683 f->bc_coeffs->ac,
2684 f->bc_coeffs->bc};
2685
2686 for (c_id = 0; c_id < 8; c_id++) {
2687 if (p[c_id] != NULL) {
2688 coeff_p[c_id] = 1;
2689 /* avoid double reads/writes in case of aliasing */
2690 for (int i = 0; i < c_id; i++) {
2691 if (p[i] == p[c_id])
2692 coeff_p[c_id] = 0;
2693 }
2694 }
2695 }
2696
2697 cs_parall_max(8, CS_INT32, coeff_p);
2698
2699 if (f->dim > 1 && coupled_key_id > -1)
2700 coupled = cs_field_get_key_int(f, coupled_key_id);
2701
2702 for (c_id = 0; c_id < 8; c_id++) {
2703
2704 char *sec_name = NULL;
2705
2706 cs_real_t *c = p[c_id];
2707
2708 if (coeff_p[c_id] == 0)
2709 continue;
2710
2711 if (coupled) {
2712 if (c_id %2 == 0)
2713 n_loc_vals = f->dim;
2714 else
2715 n_loc_vals = f->dim * f->dim;
2716 }
2717 else { /* uncoupled */
2718 n_loc_vals = f->dim;
2719 }
2720
2721 BFT_MALLOC(sec_name,
2722 strlen(f->name) + strlen(_coeff_name[c_id]) + 3,
2723 char);
2724 sprintf(sec_name, "%s::%s", f->name, _coeff_name[c_id]);
2725
2726 cs_restart_write_section(r,
2727 sec_name,
2728 3, /* location_id */
2729 n_loc_vals,
2730 CS_TYPE_cs_real_t,
2731 c);
2732
2733 BFT_FREE(sec_name);
2734
2735 if (c != p[c_id])
2736 BFT_FREE(c);
2737
2738 } /* End of loop in i (coeff type) */
2739
2740 } /* End for field with BC coeffs */
2741
2742 } /* End of loop on fields */
2743
2744 bft_printf(_(" Wrote boundary condition coefficients to checkpoint: %s\n"),
2745 cs_restart_get_name(r));
2746 }
2747
2748 /*----------------------------------------------------------------------------*/
2749 /*!
2750 * \brief Read field values from checkpoint.
2751 *
2752 * If the values are not found using the default rules based on the
2753 * field's name, its name itself, or a "restart_rename" keyed string value,
2754 * an old name may be used for compatibility with older files.
2755 * For cell-based fields, the old name base is appended automatically with
2756 * "_ce_phase01", except for scalars, where the name uses a different scheme,
2757 * based on "scalaire_ce_%04" % s_num;
2758 *
2759 * \param[in, out] r associated restart file pointer
2760 * \param[in] f_id field id
2761 * \param[in] t_id time id (0 for current, 1 for previous, ...)
2762 *
2763 * \return CS_RESTART_SUCCESS in case of success, CS_RESTART_ERR_... otherwise
2764 */
2765 /*----------------------------------------------------------------------------*/
2766
2767 int
cs_restart_read_field_vals(cs_restart_t * r,int f_id,int t_id)2768 cs_restart_read_field_vals(cs_restart_t *r,
2769 int f_id,
2770 int t_id)
2771 {
2772 cs_field_t *f = cs_field_by_id(f_id);
2773 char sec_name[128], ref_sec_name[128];
2774
2775 const char *r_name = NULL;
2776
2777 int retcode = CS_RESTART_SUCCESS;
2778
2779 /* Check for renaming */
2780
2781 int kr = cs_field_key_id_try("restart_name");
2782 if (kr > -1)
2783 r_name = cs_field_get_key_str(f, kr);
2784
2785 if (r_name == NULL)
2786 r_name = f->name;
2787
2788 /* Check for data; data will be read later, so that compatibility
2789 checks may be done first; we really try reading the data only
2790 at the end, so if it is not found, a warning will be logged only
2791 once (and not once per test), and preferentially use the
2792 base (non-compatibility) name. */
2793
2794 if (cs_restart_is_from_ncfd()) {
2795 if (strcmp(r_name, "velocity") == 0)
2796 snprintf(sec_name, 127, "%s_1", r_name);
2797 else if (strcmp(r_name, "enthalpy") == 0)
2798 snprintf(sec_name, 127, "%s_1", r_name);
2799 else if (strcmp(r_name, "temperature") == 0)
2800 snprintf(sec_name, 127, "%s_1", r_name);
2801 else if (strcmp(r_name, "k") == 0)
2802 snprintf(sec_name, 127, "%s_1", r_name);
2803 else if (strcmp(r_name, "epsilon") == 0)
2804 snprintf(sec_name, 127, "%s_1", r_name);
2805 else if (strcmp(r_name, "rij") == 0)
2806 snprintf(sec_name, 127, "%s_1", r_name);
2807 else if (strcmp(r_name, "pressure") == 0)
2808 snprintf(sec_name, 127, "%s", r_name);
2809 else
2810 snprintf(sec_name, 127, "%s::vals::%d", r_name, t_id);
2811
2812 }
2813 else
2814 snprintf(sec_name, 127, "%s::vals::%d", r_name, t_id);
2815
2816 sec_name[127] = '\0';
2817 strncpy(ref_sec_name, sec_name, 128);
2818
2819 retcode = cs_restart_check_section(r,
2820 sec_name,
2821 f->location_id,
2822 f->dim,
2823 CS_TYPE_cs_real_t);
2824
2825 /* Otherwise, try reading with basic (restart) name only if requested */
2826
2827 if ( (retcode == CS_RESTART_ERR_EXISTS || retcode == CS_RESTART_ERR_N_VALS)
2828 && r_name != f->name) {
2829 snprintf(sec_name, 127, "%s", r_name);
2830 sec_name[127] = '\0';
2831 retcode = cs_restart_check_section(r,
2832 sec_name,
2833 f->location_id,
2834 f->dim,
2835 CS_TYPE_cs_real_t);
2836
2837 } else if ((retcode == CS_RESTART_ERR_EXISTS || retcode == CS_RESTART_ERR_N_VALS)
2838 && cs_restart_is_from_ncfd() ) {
2839 /* If restart from NCFD, ensure that we can read scalar fields! */
2840 snprintf(sec_name, 127, "%s", r_name);
2841 sec_name[127] = '\0';
2842 retcode = cs_restart_check_section(r,
2843 sec_name,
2844 f->location_id,
2845 f->dim,
2846 CS_TYPE_cs_real_t);
2847 }
2848
2849 /* Read data if found */
2850
2851 if ((cs_restart_is_from_ncfd() && strcmp(f->name, "pressure") != 0) ||
2852 cs_restart_is_from_ncfd() == 0 ) {
2853
2854 if (retcode == CS_RESTART_SUCCESS)
2855 retcode = cs_restart_read_section(r,
2856 sec_name,
2857 f->location_id,
2858 f->dim,
2859 CS_TYPE_cs_real_t,
2860 f->vals[t_id]);
2861
2862 /* Try reading in compatibility mode if not found */
2863
2864 else if ( retcode == CS_RESTART_ERR_EXISTS
2865 || retcode == CS_RESTART_ERR_N_VALS) {
2866
2867 retcode = _read_field_vals_legacy(r, r_name, t_id, f);
2868
2869 /* if data is still not found, try reading in normal mode,
2870 so as to have warning/error messages relative to the
2871 current naming scheme. */
2872
2873 if ( retcode == CS_RESTART_ERR_EXISTS
2874 || retcode == CS_RESTART_ERR_N_VALS)
2875 retcode = cs_restart_read_section(r,
2876 ref_sec_name,
2877 f->location_id,
2878 f->dim,
2879 CS_TYPE_cs_real_t,
2880 f->vals[t_id]);
2881
2882 }
2883
2884 }
2885 else {
2886 /* Special case for pressure:
2887 * In NCFD it is the total pressure which is stored. Hence we need
2888 * to compute P = P_tot - Phydro - (P0-Pred0)
2889 */
2890 retcode = cs_restart_read_section(r,
2891 sec_name,
2892 f->location_id,
2893 f->dim,
2894 CS_TYPE_cs_real_t,
2895 f->vals[t_id]);
2896
2897 snprintf(sec_name, 127, "%s", "hydro_pressure");
2898 sec_name[127] = '\0';
2899 const cs_lnum_t n_cells = cs_glob_mesh->n_cells;
2900 cs_real_t *v_tmp;
2901 BFT_MALLOC(v_tmp, n_cells, cs_real_t);
2902
2903 retcode = cs_restart_read_section(r,
2904 sec_name,
2905 f->location_id,
2906 f->dim,
2907 CS_TYPE_cs_real_t,
2908 v_tmp);
2909
2910 for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
2911 f->vals[t_id][c_id] -= v_tmp[c_id]
2912 + cs_glob_fluid_properties->p0
2913 - cs_glob_fluid_properties->pred0;
2914 }
2915
2916 BFT_FREE(v_tmp);
2917 }
2918
2919 return retcode;
2920 }
2921
2922 /*----------------------------------------------------------------------------*/
2923 /*!
2924 * \brief Write field values to checkpoint.
2925 *
2926 * \param[in, out] r associated restart file pointer
2927 * \param[in] f_id field id
2928 * \param[in] t_id time id (0 for current, 1 for previous, ...)
2929 */
2930 /*----------------------------------------------------------------------------*/
2931
2932 void
cs_restart_write_field_vals(cs_restart_t * r,int f_id,int t_id)2933 cs_restart_write_field_vals(cs_restart_t *r,
2934 int f_id,
2935 int t_id)
2936 {
2937 cs_field_t *f = cs_field_by_id(f_id);
2938 char sec_name[128];
2939
2940 snprintf(sec_name, 127, "%s::vals::%d", f->name, t_id);
2941
2942 cs_restart_write_section(r,
2943 sec_name,
2944 f->location_id,
2945 f->dim,
2946 CS_TYPE_cs_real_t,
2947 f->vals[t_id]);
2948 }
2949
2950 /*----------------------------------------------------------------------------*/
2951 /*!
2952 * \brief Read restart time step info.
2953 *
2954 * \param[in, out] r associated restart file pointer
2955 */
2956 /*----------------------------------------------------------------------------*/
2957
2958 void
cs_restart_read_time_step_info(cs_restart_t * r)2959 cs_restart_read_time_step_info(cs_restart_t *r)
2960 {
2961 int retval;
2962 int _n_ts = -1;
2963 cs_real_t _ts = -1;
2964
2965 /* First syntax */
2966
2967 retval = cs_restart_check_section(r,
2968 "nbre_pas_de_temps",
2969 0,
2970 1,
2971 CS_TYPE_int);
2972
2973 if (retval == CS_RESTART_SUCCESS) {
2974 retval = cs_restart_read_section(r,
2975 "nbre_pas_de_temps",
2976 0,
2977 1,
2978 CS_TYPE_int,
2979 &_n_ts);
2980 if (retval == CS_RESTART_SUCCESS)
2981 retval = cs_restart_read_section(r,
2982 "instant_precedent",
2983 0,
2984 1,
2985 CS_TYPE_cs_real_t,
2986 &_ts);
2987
2988 if (retval == CS_RESTART_SUCCESS)
2989 cs_time_step_define_prev(_n_ts, _ts);
2990
2991 return;
2992 }
2993
2994 /* Second syntax */
2995
2996 retval = cs_restart_check_section(r,
2997 "ntcabs",
2998 0,
2999 1,
3000 CS_TYPE_int);
3001
3002 if (retval == CS_RESTART_SUCCESS) {
3003 retval = cs_restart_read_section(r,
3004 "ntcabs",
3005 0,
3006 1,
3007 CS_TYPE_int,
3008 &_n_ts);
3009 if (retval == CS_RESTART_SUCCESS)
3010 retval = cs_restart_read_section(r,
3011 "ttcabs",
3012 0,
3013 1,
3014 CS_TYPE_cs_real_t,
3015 &_ts);
3016
3017 if (retval == CS_RESTART_SUCCESS)
3018 cs_time_step_define_prev(_n_ts, _ts);
3019 return;
3020 }
3021 }
3022
3023 /*----------------------------------------------------------------------------*/
3024 /*!
3025 * \brief Loop over all fields and save them in the restart file which id is
3026 * passed in argument if it matches their "restart_file" key value.
3027 *
3028 * \param[in, out] r associated restart file pointer
3029 * \param[in] r_id value of the key "restart_file"
3030 */
3031 /*----------------------------------------------------------------------------*/
3032
3033 void
cs_restart_write_fields(cs_restart_t * r,cs_restart_file_t r_id)3034 cs_restart_write_fields(cs_restart_t *r,
3035 cs_restart_file_t r_id)
3036 {
3037 int n_fields = cs_field_n_fields();
3038 const int restart_file_key_id = cs_field_key_id("restart_file");
3039
3040 for (int f_id = 0; f_id < n_fields; f_id++) {
3041 cs_field_t *f = cs_field_by_id(f_id);
3042 if (cs_field_get_key_int(f, restart_file_key_id) == r_id) {
3043 cs_restart_write_field_vals(r, f_id, 0);
3044 }
3045 }
3046 }
3047
3048 /*----------------------------------------------------------------------------*/
3049 /*!
3050 * \brief Loop over all fields and read them in the restart file which id is
3051 * passed in argument if it matches their "restart_file" key value.
3052 *
3053 * \param[in, out] r associated restart file pointer
3054 * \param[in] r_id value of the key "restart_file"
3055 */
3056 /*----------------------------------------------------------------------------*/
3057
3058 void
cs_restart_read_fields(cs_restart_t * r,cs_restart_file_t r_id)3059 cs_restart_read_fields(cs_restart_t *r,
3060 cs_restart_file_t r_id)
3061 {
3062 int n_fields = cs_field_n_fields();
3063 const int restart_file_key_id = cs_field_key_id("restart_file");
3064
3065 for (int f_id = 0; f_id < n_fields; f_id++) {
3066 cs_field_t *f = cs_field_by_id(f_id);
3067 if (cs_field_get_key_int(f, restart_file_key_id) == r_id) {
3068 cs_restart_read_field_vals(r, f_id, 0);
3069 }
3070 }
3071 }
3072
3073 /*----------------------------------------------------------------------------*/
3074
3075 END_C_DECLS
3076