1 /*============================================================================
2 * Base time step data.
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 <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_error.h"
46 #include "bft_printf.h"
47
48 #include "cs_base.h"
49 #include "cs_log.h"
50 #include "cs_map.h"
51 #include "cs_parall.h"
52
53 /*----------------------------------------------------------------------------
54 * Header for the current file
55 *----------------------------------------------------------------------------*/
56
57 #include "cs_time_step.h"
58
59 /*----------------------------------------------------------------------------*/
60
61 BEGIN_C_DECLS
62
63 /*=============================================================================
64 * Additional doxygen documentation
65 *============================================================================*/
66
67 /*!
68 \file cs_time_step.c
69 base time step data.
70
71 \enum cs_time_step_type_t
72
73 \brief Time step type
74
75 \var CS_TIME_STEP_STEADY
76 Steady (SIMPLE) algorithm
77 \var CS_TIME_STEP_CONSTANT
78 Unsteady (SIMPLEC) algorithm with constant time step
79 \var CS_TIME_STEP_ADAPTIVE
80 Unsteady (SIMPLEC) algorithm with time-adaptive (varying) time step
81 \var CS_TIME_STEP_LOCAL
82 Pseudo-steady with time-and-space varying time step (SIMPLEC)
83
84 \struct cs_time_step_t
85
86 \brief time step descriptor
87
88 Members of this time step are publicly accessible, to allow for concise
89 syntax, as it is expected to be used in many places.
90
91 \var cs_time_step_t::is_variable
92 0 if time step is fixed in time, 1 otherwise
93 \var cs_time_step_t::is_local
94 0 if time step is uniform in space, 0 if is is local
95 (in which case the time value is only a reference)
96 \var cs_time_step_t::nt_prev
97 Absolute time step number for previous calculation.
98 In the case of a restart calculation, \ref nt_prev
99 is read from the restart file. Otherwise, it is
100 initialised to 0 \ref nt_prev is initialised
101 automatically by the code, its value is not to be
102 modified by the user.
103 \var cs_time_step_t::nt_cur
104 current absolute time step number
105 In case of restart, this is equal to \ref nt_prev + number
106 of new iterations.
107 \var cs_time_step_t::nt_max
108 maximum absolute time step number
109 For the restart calculations, \ref nt_max takes into
110 account the number of time steps of the previous calculations.
111 For instance, after a first calculation of 3 time steps, a
112 restart file of 2 time steps is realised by setting
113 \ref nt_max = 3+2 = 5
114 \var cs_time_step_t::nt_ini
115 number of time step for initialization
116 \var cs_time_step_t::t_prev
117 Absolute time value for previous calculation.
118 In the case of a restart calculation, \ref t_prev is read from
119 the restart file. Otherwise it is initialised to 0.\n
120 \ref t_prev is initialised automatically by the code,
121 its value is not to be modified by the user.
122 \var cs_time_step_t::t_cur
123 Current absolute time.
124 For the restart calculations, \ref t_cur takes
125 into account the physical time of the previous calculations.\n
126 If the time step is uniform (\ref cs_time_step_options_t::idtvar "idtvar"
127 = CS_TIME_STEP_CONSTANT or CS_ADAPTATIVE_TIME_STEP),
128 \ref t_cur increases of \ref dt (value of the time step)
129 at each iteration.
130 If the time step is non-uniform (\ref cs_time_step_options_t::idtvar
131 "idtvar"=CS_TIME_STEP_LOCAL),
132 \ref t_cur increases of \ref cs_time_step_t::dt_ref "dt_ref" at each
133 time step.\n
134 \ref t_cur is initialised and updated automatically by the code,
135 its value is not to be modified by the user.
136 \var cs_time_step_t::t_max
137 maximum absolute time
138
139 \struct cs_time_step_options_t
140
141 \brief time step options descriptor
142
143 Members of this time step options descriptor are publicly accessible, to
144 allow for concise syntax.
145
146 \var cs_time_step_options_t::iptlro
147 Clip the time step with respect to the buoyant effects\n
148
149 When density gradients and gravity are present, a local thermal time
150 step can be calculated, based on the Brunt-Vaisala frequency. In
151 numerical simulations, it is usually wise for the time step to be
152 lower than this limit, otherwise numerical instabilities may appear.\n
153 \ref iptlro indicates whether the time step should be limited to the
154 local thermal time step (=1) or not (=0).\n
155 When \ref iptlro=1, the log shows the number of cells where the
156 time step has been clipped due to the thermal criterion, as well as
157 the maximum ratio between the time step and the maximum thermal time
158 step. If \ref idtvar=CS_TIME_STEP_CONSTANT, since the time step is fixed
159 and cannot be clipped, this ratio can be greater than 1.
160 When \ref idtvar > CS_TIME_STEP_CONSTANT, this
161 ratio will be less than 1, except if the constraint \ref dtmin has
162 prevented the code from reaching a sufficiently low value for \ref dt.
163 Useful when density gradients and gravity are present.
164
165 \var cs_time_step_options_t::idtvar
166 Time step type: constant, adaptive, or steady algorithm variants.
167 If the numerical scheme is a second-order in time, only the
168 option CS_TIME_STEP_CONSTANT is allowed.
169 \anchor idtvar
170
171 \var cs_time_step_t::dt_ref
172 Reference time step.\n
173
174 This is the time step value used in the case of a calculation run with a
175 uniform and constant time step, i.e. \ref idtvar =CS_TIME_STEP_CONSTANT
176 (restart calculation or not). It is the value used to initialize the
177 time step in the case of an initial calculation run with a non-constant
178 time step (\ref idtvar=CS_TIME_STEP_ADAPTATIVE or CS_TIME_STEP_LOCAL).
179 It is also the value used to initialise the time step in the case of
180 a restart calculation in which the type of time step has been changed
181 (for instance, \ref idtvar=CS_TIME_STEP_ADAPTATIVE in the new calculation
182 and \ref idtvar = CS_TIME_STEP_CONSTANT or CS_STEP_LOCAL_TIME in the
183 previous calculation).\n
184 See \ref user_initialization_time_step for examples.
185
186 \var cs_time_step_options_t::coumax
187 Maximum Courant number (when \ref idtvar is different
188 from CS_TIME_STEP_CONSTANT).
189
190 \var cs_time_step_options_t::cflmmx
191 Max. Courant number for the continuity equation in compressible model.
192
193 \var cs_time_step_options_t::foumax
194 Maximum Fourier number (when \ref idtvar is different from
195 CS_TIME_STEP_CONSTANT).
196
197 \var cs_time_step_options_t::varrdt
198 Maximum allowed relative increase in the calculated time step value
199 between two successive time steps (to ensure stability, any decrease
200 in the time step is immediate and without limit).\n
201 Useful when idtvar is different from CS_TIME_STEP_CONSTANT.
202
203 \var cs_time_step_options_t::dtmin
204 Lower limit for the calculated time step when idtvar is different from
205 CS_TIME_STEP_CONSTANT.
206 Take \ref dtmin = min (ld/ud, sqrt(lt/(gdelta rho/rho)), ...).
207
208 \var cs_time_step_options_t::dtmax
209 Upper limit for the calculated time step when idtvar is different from
210 CS_TIME_STEP_CONSTANT.
211 Take \ref dtmax = max (ld/ud, sqrt(lt/(gdelta rho/rho)), ...).
212
213 \var cs_time_step_options_t::relxst
214 Relaxation coefficient for the steady algorithm.
215 \ref relxst = 1 : no relaxation.
216
217 */
218
219 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
220
221 /*=============================================================================
222 * Macro definitions
223 *============================================================================*/
224
225 /*============================================================================
226 * Type definitions
227 *============================================================================*/
228
229 /*============================================================================
230 * Static global variables
231 *============================================================================*/
232
233 static const char *cs_time_step_type_enum_name[]
234 = {"CS_TIME_STEP_STEADY",
235 "CS_TIME_STEP_CONSTANT",
236 "CS_TIME_STEP_ADAPTIVE",
237 "CS_TIME_STEP_LOCAL"};
238
239 static const char *cs_time_step_type_name[]
240 = {N_("steady algorithm"),
241 N_("constant"),
242 N_("time-adaptive"),
243 N_("local in time and space (pseudo-steady)")};
244
245 /* main time step structure and associated pointer */
246
247 static cs_time_step_t _time_step = {
248 .is_variable = 0,
249 .is_local = 0,
250 .nt_prev = 0,
251 .nt_cur = 0,
252 .nt_max = -1,
253 .nt_ini = 2,
254 .t_prev = 0.,
255 .t_cur = 0.,
256 .t_max = -1.,
257 .dt = {0.1, 0.1, 0.1},
258 .dt_ref = 0.1,
259 .dt_next = 0.1
260 };
261
262 static cs_time_step_options_t _time_step_options = {
263 .iptlro = 0,
264 .idtvar = CS_TIME_STEP_CONSTANT, /* constant time step by default */
265 .coumax = 1.,
266 .cflmmx = 0.99,
267 .foumax = 10.,
268 .varrdt = 0.1,
269 .dtmin = -1.e13,
270 .dtmax = -1.e13,
271 .relxst = 0.7 /* Not used in CDO schemes */
272 };
273
274 const cs_time_step_t *cs_glob_time_step = &_time_step;
275
276 const cs_time_step_options_t *cs_glob_time_step_options = &_time_step_options;
277
278 static double _c = 0; /* compensation term for Kahan sum */
279
280 /*============================================================================
281 * Prototypes for functions intended for use only by Fortran wrappers.
282 * (descriptions follow, with function bodies).
283 *============================================================================*/
284
285 void
286 cs_f_time_step_get_pointers(int **nt_prev,
287 int **nt_cur,
288 int **nt_max,
289 int **nt_ini,
290 double **dtref,
291 double **t_prev,
292 double **t_cur,
293 double **t_max);
294
295 void
296 cs_f_time_step_options_get_pointers(int **iptlro,
297 int **idtvar,
298 double **coumax,
299 double **cflmmx,
300 double **foumax,
301 double **varrdt,
302 double **dtmin,
303 double **dtmax,
304 double **relxst);
305
306 /*============================================================================
307 * Private function definitions
308 *============================================================================*/
309
310 /*============================================================================
311 * Fortran wrapper function definitions
312 *============================================================================*/
313
314 /*----------------------------------------------------------------------------
315 * Get pointers to members of the global time step structure.
316 *
317 * This function is intended for use by Fortran wrappers, and
318 * enables mapping to Fortran global pointers.
319 *
320 * parameters:
321 * nt_prev --> pointer to cs_glob_time_step->nt_prev
322 * nt_cur --> pointer to cs_glob_time_step->nt_cur
323 * nt_max --> pointer to cs_glob_time_step->nt_max
324 * nt_ini --> pointer to cs_glob_time_step->nt_ini
325 * dt_ref --> pointer to cs_glob_time_step->dt_ref
326 * t_prev --> pointer to cs_glob_time_step->t_prev
327 * t_cur --> pointer to cs_glob_time_step->t_cur
328 * t_max --> pointer to cs_glob_time_step->t_ax
329 *----------------------------------------------------------------------------*/
330
331 void
cs_f_time_step_get_pointers(int ** nt_prev,int ** nt_cur,int ** nt_max,int ** nt_ini,double ** dtref,double ** t_prev,double ** t_cur,double ** t_max)332 cs_f_time_step_get_pointers(int **nt_prev,
333 int **nt_cur,
334 int **nt_max,
335 int **nt_ini,
336 double **dtref,
337 double **t_prev,
338 double **t_cur,
339 double **t_max)
340 {
341 *nt_prev = &(_time_step.nt_prev);
342 *nt_cur = &(_time_step.nt_cur);
343 *nt_max = &(_time_step.nt_max);
344 *nt_ini = &(_time_step.nt_ini);
345 *dtref = &(_time_step.dt_ref);
346 *t_prev = &(_time_step.t_prev);
347 *t_cur = &(_time_step.t_cur);
348 *t_max = &(_time_step.t_max);
349 }
350
351 /*----------------------------------------------------------------------------
352 * Get pointers to members of the global time step structure.
353 *
354 * This function is intended for use by Fortran wrappers, and
355 * enables mapping to Fortran global pointers.
356 *
357 * parameters:
358 * iptlro --> pointer to cs_glob_time_step_options->iptlro
359 * idtvar --> pointer to cs_glob_time_step_options->idtvar
360 * coumax --> pointer to cs_glob_time_step_options->coumax
361 * cflmmx --> pointer to cs_glob_time_step_options->cflmmx
362 * foumax --> pointer to cs_glob_time_step_options->foumax
363 * varrdt --> pointer to cs_glob_time_step_options->varrdt
364 * dtmin --> pointer to cs_glob_time_step_options->dtmin
365 * dtmax --> pointer to cs_glob_time_step_options->dtmax
366 * relxst --> pointer to cs_glob_time_step_options->relxst
367 *----------------------------------------------------------------------------*/
368
369 void
cs_f_time_step_options_get_pointers(int ** iptlro,int ** idtvar,double ** coumax,double ** cflmmx,double ** foumax,double ** varrdt,double ** dtmin,double ** dtmax,double ** relxst)370 cs_f_time_step_options_get_pointers(int **iptlro,
371 int **idtvar,
372 double **coumax,
373 double **cflmmx,
374 double **foumax,
375 double **varrdt,
376 double **dtmin,
377 double **dtmax,
378 double **relxst)
379 {
380 *iptlro = &(_time_step_options.iptlro);
381 *idtvar = &(_time_step_options.idtvar);
382 *coumax = &(_time_step_options.coumax);
383 *cflmmx = &(_time_step_options.cflmmx);
384 *foumax = &(_time_step_options.foumax);
385 *varrdt = &(_time_step_options.varrdt);
386 *dtmin = &(_time_step_options.dtmin );
387 *dtmax = &(_time_step_options.dtmax );
388 *relxst = &(_time_step_options.relxst);
389 }
390
391 /*=============================================================================
392 * Private function definitions
393 *============================================================================*/
394
395 /*----------------------------------------------------------------------------*/
396 /*!
397 * \brief Update Kahan compensation.
398 *
399 * Called as a function to help avoid compiler optimizating Kahan
400 * summation out (though aggressive or LTO optimizations might).
401 *
402 * new copensation term is (a - b) - c
403 *
404 * \param[in] a a
405 * \param[in] b b
406 * \param[in] c c
407 */
408 /*----------------------------------------------------------------------------*/
409
410 static void
_update_kahan_compensation(double * a,double * b,double * c)411 _update_kahan_compensation(double *a,
412 double *b,
413 double *c)
414 {
415 _c = (*a - *b) - *c;
416 }
417
418 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
419
420 /*=============================================================================
421 * Public function definitions
422 *============================================================================*/
423
424 /*----------------------------------------------------------------------------*/
425 /*!
426 * \brief Provide read/write access to cs_glob_time_step
427 *
428 * \return pointer to global time step structure
429 */
430 /*----------------------------------------------------------------------------*/
431
432 cs_time_step_t *
cs_get_glob_time_step(void)433 cs_get_glob_time_step(void)
434 {
435 return &_time_step;
436 }
437
438 /*----------------------------------------------------------------------------*/
439 /*!
440 * \brief Provide read/write access to cs_glob_time_step_options
441 *
442 * \return pointer to global time step options structure
443 */
444 /*----------------------------------------------------------------------------*/
445
446 cs_time_step_options_t *
cs_get_glob_time_step_options(void)447 cs_get_glob_time_step_options(void)
448 {
449 return &_time_step_options;
450 }
451
452 /*----------------------------------------------------------------------------*/
453 /*!
454 * \brief Define whether time step is variable or not
455 *
456 * \param[in] is_variable 0 if time step is variable in time, 1 if it is fixed
457 */
458 /*----------------------------------------------------------------------------*/
459
460 void
cs_time_step_define_variable(int is_variable)461 cs_time_step_define_variable(int is_variable)
462 {
463 if (is_variable > 0)
464 _time_step.is_variable = 1;
465 else
466 _time_step.is_variable = 0;
467 }
468
469 /*----------------------------------------------------------------------------*/
470 /*!
471 * \brief Define whether time step is local in space or not
472 *
473 * \param[in] is_local 0 if time step is uniform in space, 1 if it is local
474 */
475 /*----------------------------------------------------------------------------*/
476
477 void
cs_time_step_define_local(int is_local)478 cs_time_step_define_local(int is_local)
479 {
480 if (is_local > 0)
481 _time_step.is_local = 1;
482 else
483 _time_step.is_local = 0;
484 }
485
486 /*----------------------------------------------------------------------------*/
487 /*!
488 * \brief Define maximum time step number
489 *
490 * \param[in] nt_max maximum time step number (unlimited if negative)
491 */
492 /*----------------------------------------------------------------------------*/
493
494 void
cs_time_step_define_nt_max(int nt_max)495 cs_time_step_define_nt_max(int nt_max)
496 {
497 _time_step.nt_max = nt_max;
498 _time_step.t_max = -1.;
499 }
500
501 /*----------------------------------------------------------------------------*/
502 /*!
503 * \brief Define maximum time value
504 *
505 * \param[in] t_max maximum time value (unlimited if negative)
506 */
507 /*----------------------------------------------------------------------------*/
508
509 void
cs_time_step_define_t_max(double t_max)510 cs_time_step_define_t_max(double t_max)
511 {
512 _time_step.t_max = t_max;
513 _time_step.nt_max = -1;
514 }
515
516 /*----------------------------------------------------------------------------*/
517 /*!
518 * \brief Set time values from previous (usually restarted) calculations
519 *
520 * \param[in] nt_prev previous time step number
521 * \param[in] t_prev previous physical time
522 */
523 /*----------------------------------------------------------------------------*/
524
525 void
cs_time_step_define_prev(int nt_prev,double t_prev)526 cs_time_step_define_prev(int nt_prev,
527 double t_prev)
528 {
529 _time_step.nt_prev = nt_prev;
530 _time_step.nt_cur = nt_prev;
531 _time_step.t_prev = t_prev;
532 _time_step.t_cur = t_prev;
533 }
534
535 /*----------------------------------------------------------------------------*/
536 /*!
537 * \brief Increment the global time step.
538 *
539 * \param[in] dt time step value to increment
540 */
541 /*----------------------------------------------------------------------------*/
542
543 void
cs_time_step_increment(double dt)544 cs_time_step_increment(double dt)
545 {
546 _time_step.dt[2] = _time_step.dt[1];
547 _time_step.dt[1] = _time_step.dt[0];
548 _time_step.dt[0] = dt;
549
550 double z = dt - _c;
551 double t = _time_step.t_cur + z;
552
553 _update_kahan_compensation(&t, &_time_step.t_cur, &z);
554
555 _time_step.t_cur = t;
556 _time_step.nt_cur += 1;
557
558 if (_time_step.is_local)
559 cs_base_update_status("time step: %d\n",
560 _time_step.nt_cur);
561 else
562 cs_base_update_status("time step: %d; t = %g\n",
563 _time_step.nt_cur, _time_step.t_cur);
564 }
565
566 /*----------------------------------------------------------------------------*/
567 /*!
568 * \brief Redefine the current time values.
569 *
570 * \remark Using \ref cs_time_step_increment is preferred, but this function
571 * may be required for reverting to a previous time step.
572 *
573 * \param[in] nt_cur current time step number
574 * \param[in] t_cur current physical time
575 */
576 /*----------------------------------------------------------------------------*/
577
578 void
cs_time_step_redefine_cur(int nt_cur,double t_cur)579 cs_time_step_redefine_cur(int nt_cur,
580 double t_cur)
581 {
582 _time_step.nt_cur = nt_cur;
583 _time_step.t_cur = t_cur;
584 }
585
586 /*----------------------------------------------------------------------------*/
587 /*!
588 * \brief Print the time stepping options to setup.log.
589 */
590 /*----------------------------------------------------------------------------*/
591
592 void
cs_time_step_log_setup(void)593 cs_time_step_log_setup(void)
594 {
595 cs_log_printf
596 (CS_LOG_SETUP,
597 _("\nTime stepping options\n"
598 "---------------------\n\n"));
599
600 cs_time_step_type_t ts_type = cs_glob_time_step_options->idtvar;
601 int ts_id = ts_type + 1; /* (this enum starts at -1) */
602
603 if (ts_type == CS_TIME_STEP_STEADY) {
604
605 /* Relaxation parameters */
606 cs_log_printf
607 (CS_LOG_SETUP,
608 _(" Steady (SIMPLE) algorithm\n\n"
609 " Global parameters\n\n"
610 " idtvar: %21s (%s)\n"
611 " relxst: %17.5g (Reference relaxation coefficient)\n\n"),
612 cs_time_step_type_enum_name[ts_id],
613 cs_time_step_type_name[ts_id],
614 cs_glob_time_step_options->relxst);
615 }
616
617 else {
618
619 if (ts_type ==CS_TIME_STEP_CONSTANT) {
620 cs_log_printf
621 (CS_LOG_SETUP,
622 _(" Unsteady algorithm\n\n"
623 " Time step parameters\n\n"
624 " idtvar: %21s (%s)\n"
625 " dtref: %17.5g (Reference time step)\n\n"),
626 cs_time_step_type_enum_name[ts_id],
627 cs_time_step_type_name[ts_id],
628 cs_glob_time_step->dt_ref);
629 }
630
631 else {
632 if (cs_glob_time_step_options->idtvar == CS_TIME_STEP_ADAPTIVE) {
633 cs_log_printf
634 (CS_LOG_SETUP,
635 _(" Unsteady algorithm\n\n"));
636 } else if (cs_glob_time_step_options->idtvar == CS_TIME_STEP_LOCAL) {
637 cs_log_printf
638 (CS_LOG_SETUP,
639 _(" Space & time varying time step algorithm (pseudo-steady)\n\n"));
640 }
641
642 cs_log_printf
643 (CS_LOG_SETUP,
644 _(" Time step parameters:\n\n"
645 " idtvar: %21s (%s)\n"
646 " iptlro: %17d (1: rho-related DT clipping)\n"
647 " coumax: %17.5g (Maximum target CFL)\n"
648 " foumax: %17.5g (Maximum target Fourier)\n"
649 " varrdt: %17.5g (For var. DT, max. increase)\n"
650 " dtmin: %17.5g (Minimum time step)\n"
651 " dtmax: %17.5g (Maximum time step)\n"
652 " dtref: %17.5g (Reference time step)\n\n"
653 " When the value of coumax or foumax is negative\n"
654 " or zero, the associated time step limitation\n"
655 " (for CFL and Fourier respectively) is ignored.\n\n"),
656 cs_time_step_type_enum_name[ts_id],
657 cs_time_step_type_name[ts_id],
658 cs_glob_time_step_options->iptlro,
659 cs_glob_time_step_options->coumax,
660 cs_glob_time_step_options->foumax,
661 cs_glob_time_step_options->varrdt,
662 cs_glob_time_step_options->dtmin,
663 cs_glob_time_step_options->dtmax,
664 cs_glob_time_step->dt_ref);
665
666 }
667 }
668 }
669
670 /*----------------------------------------------------------------------------*/
671
672 END_C_DECLS
673