ryujin 2.1.1 revision 94ce4289f7944859085a4b6724bdcb05a682328b
time_loop.template.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
3// Copyright (C) 2020 - 2024 by the ryujin authors
4//
5
6#pragma once
7
8#include "instrumentation.h"
9#include "scope.h"
10#include "state_vector.h"
11#include "time_loop.h"
12#include "version_info.h"
13
14#include <deal.II/base/logstream.h>
15#include <deal.II/base/work_stream.h>
16#include <deal.II/numerics/vector_tools.h>
17#include <deal.II/numerics/vector_tools.templates.h>
18
19#include <filesystem>
20#include <fstream>
21#include <iomanip>
22
23using namespace dealii;
24
25namespace ryujin
26{
27 template <typename Description, int dim, typename Number>
29 : ParameterAcceptor("/A - TimeLoop")
30 , mpi_ensemble_(mpi_comm)
31 , hyperbolic_system_(mpi_ensemble_, "/B - Equation")
32 , parabolic_system_(mpi_ensemble_, "/B - Equation")
33 , discretization_(mpi_ensemble_, "/C - Discretization")
34 , offline_data_(mpi_ensemble_, discretization_, "/D - OfflineData")
35 , initial_values_(mpi_ensemble_,
36 "/E - InitialValues",
37 mpi_ensemble_,
38 offline_data_,
39 hyperbolic_system_,
40 parabolic_system_)
41 , hyperbolic_module_(mpi_ensemble_,
42 computing_timer_,
43 offline_data_,
44 hyperbolic_system_,
45 initial_values_,
46 "/F - HyperbolicModule")
47 , parabolic_module_(mpi_ensemble_,
48 computing_timer_,
49 offline_data_,
50 hyperbolic_system_,
51 parabolic_system_,
52 initial_values_,
53 "/G - ParabolicModule")
54 , time_integrator_(mpi_ensemble_,
55 offline_data_,
56 hyperbolic_module_,
57 parabolic_module_,
58 "/H - TimeIntegrator")
59 , mesh_adaptor_(mpi_ensemble_,
60 offline_data_,
61 hyperbolic_system_,
62 parabolic_system_,
63 hyperbolic_module_.initial_precomputed(),
64 hyperbolic_module_.alpha(),
65 "/I - MeshAdaptor")
66 , solution_transfer_(
67 mpi_ensemble_, offline_data_, hyperbolic_system_, parabolic_system_)
68 , postprocessor_(mpi_ensemble_,
69 offline_data_,
70 hyperbolic_system_,
71 parabolic_system_,
72 "/J - VTUOutput")
73 , vtu_output_(mpi_ensemble_,
74 offline_data_,
75 hyperbolic_system_,
76 parabolic_system_,
77 postprocessor_,
78 hyperbolic_module_.initial_precomputed(),
79 hyperbolic_module_.alpha(),
80 "/J - VTUOutput")
81 , quantities_(mpi_ensemble_,
82 offline_data_,
83 hyperbolic_system_,
84 parabolic_system_,
85 "/K - Quantities")
86 {
87 base_name_ = "test";
88 add_parameter("basename", base_name_, "Base name for all output files");
89
90 t_final_ = Number(5.);
91 add_parameter("final time", t_final_, "Final time");
92
93 enforce_t_final_ = false;
94 add_parameter("enforce final time",
95 enforce_t_final_,
96 "Boolean indicating whether the final time should be "
97 "enforced strictly. If set to true the last time step is "
98 "shortened so that the simulation ends precisely at t_final");
99
100 timer_granularity_ = Number(0.01);
101 add_parameter("timer granularity",
102 timer_granularity_,
103 "The timer granularity specifies the time interval after "
104 "which compute, output, postprocessing, and mesh adaptation "
105 "routines are run. This \"baseline tick\" is further "
106 "modified by the corresponding \"*_multiplier\" options");
107
108 enable_output_full_ = false;
109 add_parameter("enable output full",
110 enable_output_full_,
111 "Write out full pvtu records. The frequency is determined by "
112 "\"timer granularity\" and \"timer output full multiplier\"");
113
114 enable_output_levelsets_ = false;
115 add_parameter(
116 "enable output levelsets",
117 enable_output_levelsets_,
118 "Write out levelsets pvtu records. The frequency is determined by "
119 "\"timer granularity\" and \"timer output levelsets multiplier\"");
120
121 enable_compute_error_ = false;
122 add_parameter("enable compute error",
123 enable_compute_error_,
124 "Flag to control whether we compute the Linfty Linf_norm of "
125 "the difference to an analytic solution. Implemented only "
126 "for certain initial state configurations.");
127
128 enable_compute_quantities_ = false;
129 add_parameter(
130 "enable compute quantities",
131 enable_compute_quantities_,
132 "Flag to control whether we compute quantities of interest. The "
133 "frequency how often quantities are logged is determined by \"timer "
134 "granularity\" and \"timer compute quantities multiplier\"");
135
136 enable_mesh_adaptivity_ = false;
137 add_parameter(
138 "enable mesh adaptivity",
139 enable_mesh_adaptivity_,
140 "Flag to control whether we use an adaptive mesh refinement strategy. "
141 "The frequency how often we query MeshAdaptor::analyze() for deciding "
142 "on adapting the mesh is determined by \"timer granularity\" and "
143 "\"timer mesh refinement multiplier\"");
144
145 timer_output_full_multiplier_ = 1;
146 add_parameter("timer output full multiplier",
147 timer_output_full_multiplier_,
148 "Multiplicative modifier applied to \"timer granularity\" "
149 "that determines the full pvtu writeout granularity");
150
151 timer_output_levelsets_multiplier_ = 1;
152 add_parameter("timer output levelsets multiplier",
153 timer_output_levelsets_multiplier_,
154 "Multiplicative modifier applied to \"timer granularity\" "
155 "that determines the levelsets pvtu writeout granularity");
156
157 timer_compute_quantities_multiplier_ = 1;
158 add_parameter(
159 "timer compute quantities multiplier",
160 timer_compute_quantities_multiplier_,
161 "Multiplicative modifier applied to \"timer granularity\" that "
162 "determines the writeout granularity for quantities of interest");
163
164 std::copy(std::begin(View::component_names),
165 std::end(View::component_names),
166 std::back_inserter(error_quantities_));
167
168 add_parameter("error quantities",
169 error_quantities_,
170 "List of conserved quantities used in the computation of the "
171 "error norms.");
172
173 error_normalize_ = true;
174 add_parameter("error normalize",
175 error_normalize_,
176 "Flag to control whether the error should be normalized by "
177 "the corresponding norm of the analytic solution.");
178
179 resume_ = false;
180 add_parameter("resume", resume_, "Resume an interrupted computation");
181
182 resume_at_time_zero_ = false;
183 add_parameter("resume at time zero",
184 resume_at_time_zero_,
185 "Resume from the latest checkpoint but set the time to t=0.");
186
187 terminal_update_interval_ = 5;
188 add_parameter("terminal update interval",
189 terminal_update_interval_,
190 "Number of seconds after which output statistics are "
191 "recomputed and printed on the terminal. Setting the "
192 "interval to zero disables terminal output.");
193
194 terminal_show_rank_throughput_ = true;
195 add_parameter("terminal show rank throughput",
196 terminal_show_rank_throughput_,
197 "If set to true an average per rank throughput is computed "
198 "by dividing the total consumed CPU time (per rank) by the "
199 "number of threads (per rank). If set to false then a plain "
200 "average per thread \"CPU\" throughput value is computed by "
201 "using the umodified total accumulated CPU time.");
202
203 checkpoint_update_interval_ = 0;
204 add_parameter(
205 "checkpoint update interval",
206 checkpoint_update_interval_,
207 "Number of seconds after which a new checkpoint is written out to "
208 "disk. Setting the interval to zero disables checkpointing.");
209
210 debug_filename_ = "";
211 add_parameter("debug filename",
212 debug_filename_,
213 "If set to a nonempty string then we output the contents of "
214 "this file at the end. This is mainly useful in the "
215 "testsuite to output files we wish to compare");
216 }
217
218
219 /*
220 * ---------------------------------------------------------------------------
221 * Setup and main loop:
222 * ---------------------------------------------------------------------------
223 */
224
225
226 template <typename Description, int dim, typename Number>
228 {
229#ifdef DEBUG_OUTPUT
230 std::cout << "TimeLoop<dim, Number>::run()" << std::endl;
231#endif
232
233 {
234 base_name_ensemble_ = base_name_;
235 if (mpi_ensemble_.n_ensembles() > 1) {
236 print_info("setting up MPI ensemble");
237 base_name_ensemble_ += "-ensemble_" + dealii::Utilities::int_to_string(
238 mpi_ensemble_.ensemble(),
239 mpi_ensemble_.n_ensembles());
240 }
241 }
242
243 /* Attach log file and record runtime parameters: */
244
245 if (mpi_ensemble_.world_rank() == 0)
246 logfile_.open(base_name_ + ".log");
247
248 print_parameters(logfile_);
249
250 /*
251 * Prepare data structures:
252 */
253
254 Number t = 0.;
255 unsigned int timer_cycle = 0;
256 StateVector state_vector;
257
258 /* Create a small lambda for preparing compute kernels: */
259 const auto prepare_compute_kernels = [&]() {
260 print_info("preparing compute kernels");
261
262 unsigned int n_parabolic_state_vectors =
263 parabolic_system_.get().n_parabolic_state_vectors();
264
265 offline_data_.prepare(
266 problem_dimension, n_precomputed_values, n_parabolic_state_vectors);
267
268 hyperbolic_module_.prepare();
269 parabolic_module_.prepare();
270 time_integrator_.prepare();
271 mesh_adaptor_.prepare(/*needs current timepoint*/ t);
272 postprocessor_.prepare();
273 vtu_output_.prepare();
274 quantities_.prepare(base_name_ensemble_);
275 print_mpi_partition(logfile_);
276
277 if (mpi_ensemble_.ensemble_rank() == 0)
278 n_global_dofs_ = dealii::Utilities::MPI::sum(
279 offline_data_.dof_handler().n_dofs(),
280 mpi_ensemble_.ensemble_leader_communicator());
281 };
282
283 {
284 Scope scope(computing_timer_, "(re)initialize data structures");
285 print_info("initializing data structures");
286
287 if (resume_) {
288 print_info("resume: reading mesh and loading state vector");
289
290 read_checkpoint(state_vector,
291 base_name_ensemble_,
292 t,
293 timer_cycle,
294 prepare_compute_kernels);
295
296 if (resume_at_time_zero_) {
297 /* Reset the current time t and the output cycle count to zero: */
298 t = 0.;
299 timer_cycle = 0;
300 }
301
302 } else {
303 print_info("creating mesh and interpolating initial values");
304
305 discretization_.prepare(base_name_ensemble_);
306
307 prepare_compute_kernels();
308
309 Vectors::reinit_state_vector<Description>(state_vector, offline_data_);
310 std::get<0>(state_vector) =
311 initial_values_.get().interpolate_hyperbolic_vector();
312 }
313 }
314
315 /*
316 * In debug mode poison constrained degrees of freedom and precomputed
317 * values:
318 */
319 Vectors::debug_poison_constrained_dofs<Description>(state_vector,
320 offline_data_);
321 Vectors::debug_poison_precomputed_values<Description>(state_vector,
322 offline_data_);
323
324 Number last_terminal_output = terminal_update_interval_ == Number(0.)
325 ? std::numeric_limits<Number>::max()
326 : Number(0.);
327 Number last_checkpoint = checkpoint_update_interval_ == Number(0.)
328 ? std::numeric_limits<Number>::max()
329 : Number(0.);
330
331 /*
332 * The honorable main loop:
333 */
334
335 print_info("entering main loop");
336 computing_timer_["time loop"].start();
337
338 constexpr Number relax =
339 Number(1.) - Number(10.) * std::numeric_limits<Number>::epsilon();
340
341 unsigned int cycle = 1;
342 for (;; ++cycle) {
343
344#ifdef DEBUG_OUTPUT
345 std::cout << "\n\n### cycle = " << cycle << " ###\n\n" << std::endl;
346#endif
347
348 /* Accumulate quantities of interest: */
349
350 if (enable_compute_quantities_) {
351 Scope scope(computing_timer_,
352 "time step [X] - accumulate quantities");
353 quantities_.accumulate(state_vector, t);
354 }
355
356 /* Perform output tasks whenever we reach a timer tick: */
357
358 if (t >= relax * timer_cycle * timer_granularity_) {
359 if (enable_compute_error_) {
360 /*
361 * FIXME: We interpolate the analytic solution at every timer
362 * tick. If we happen to actually not output anything then this
363 * is terribly inefficient...
364 */
365
366 StateVector analytic;
367 {
368 Scope scope(computing_timer_,
369 "time step [X] - interpolate analytic solution");
370 Vectors::reinit_state_vector<Description>(analytic, offline_data_);
371 std::get<0>(analytic) =
372 initial_values_.get().interpolate_hyperbolic_vector(t);
373 }
374
375 output(analytic,
376 base_name_ensemble_ + "-analytic_solution",
377 t,
378 timer_cycle);
379 }
380
381 output(state_vector, base_name_ensemble_ + "-solution", t, timer_cycle);
382
383 if (enable_compute_quantities_ &&
384 (timer_cycle % timer_compute_quantities_multiplier_ == 0)) {
385 Scope scope(computing_timer_,
386 "time step [X] - write out quantities");
387 quantities_.write_out(state_vector, t, timer_cycle);
388 }
389
390 ++timer_cycle;
391 }
392
393 /* Break if we have reached the final time. */
394
395 if (t >= relax * t_final_)
396 break;
397
398 /* Peform a mesh adaptation cycle: */
399
400 if (enable_mesh_adaptivity_) {
401 {
402 Scope scope(computing_timer_,
403 "time step [X] - analyze for mesh adaptation");
404
405 mesh_adaptor_.analyze(state_vector, t, cycle);
406 }
407
408 if (mesh_adaptor_.need_mesh_adaptation()) {
409 Scope scope(computing_timer_, "(re)initialize data structures");
410 print_info("performing mesh adaptation");
411
412 hyperbolic_module_.prepare_state_vector(state_vector, t);
413 if (!ParabolicSystem::is_identity)
414 parabolic_module_.prepare_state_vector(state_vector, t);
415
416 adapt_mesh_and_transfer_state_vector(state_vector,
417 prepare_compute_kernels);
418 }
419 }
420
421 /* Perform a time step: */
422
423 const auto tau = time_integrator_.step(
424 state_vector,
425 t,
426 enforce_t_final_
427 ? std::min(t_final_, timer_cycle * timer_granularity_)
428 : std::numeric_limits<Number>::max());
429
430 t += tau;
431
432 /* Synchronize wall time: */
433
434 auto wall_time = computing_timer_["time loop"].wall_time();
435 {
436 Scope scope(computing_timer_,
437 "time step [X] _ - synchronization barriers");
438 wall_time =
439 Utilities::MPI::max(wall_time, mpi_ensemble_.world_communicator());
440 }
441
442 /* Print and record cycle statistics: */
443
444 const bool write_to_log_file =
445 (terminal_update_interval_ != Number(0.)) && /* suppress output */
446 (t >= relax * timer_cycle * timer_granularity_);
447
448 const bool update_terminal =
449 (wall_time >= last_terminal_output + terminal_update_interval_);
450
451 if (write_to_log_file || update_terminal) {
452 Scope scope(computing_timer_,
453 "time step [X] _ - synchronization barriers");
454 print_cycle_statistics(cycle,
455 t,
456 timer_cycle,
457 last_checkpoint,
458 /*logfile*/ write_to_log_file);
459 last_terminal_output = wall_time;
460 }
461
462 const bool update_checkpoint =
463 (wall_time >= last_checkpoint + checkpoint_update_interval_);
464
465 if (update_checkpoint) {
466 Scope scop(computing_timer_, "time step [X] - perform checkpointing");
467
468 hyperbolic_module_.prepare_state_vector(state_vector, t);
469 if (!ParabolicSystem::is_identity)
470 parabolic_module_.prepare_state_vector(state_vector, t);
471
472 print_info("scheduling checkpointing");
473 write_checkpoint(state_vector, base_name_ensemble_, t, timer_cycle);
474 last_checkpoint = wall_time;
475 }
476 } /* end of loop */
477
478 /* We have actually performed one cycle less. */
479 --cycle;
480
481 if (checkpoint_update_interval_ != Number(0.)) {
482 Scope scope(computing_timer_, "time step [X] - perform checkpointing");
483
484 hyperbolic_module_.prepare_state_vector(state_vector, t);
485 if (!ParabolicSystem::is_identity)
486 parabolic_module_.prepare_state_vector(state_vector, t);
487
488 print_info("scheduling checkpointing");
489 write_checkpoint(state_vector, base_name_ensemble_, t, timer_cycle);
490 }
491
492 computing_timer_["time loop"].stop();
493
494 if (terminal_update_interval_ != Number(0.)) {
495 /* Write final timing statistics to screen and logfile: */
496 print_cycle_statistics(cycle,
497 t,
498 timer_cycle,
499 last_checkpoint,
500 /*logfile*/ true,
501 /*final*/ true);
502 }
503
504 if (enable_compute_error_) {
505 /* Output final error: */
506 compute_error(state_vector, t);
507 }
508
509 if (mpi_ensemble_.world_rank() == 0 && debug_filename_ != "") {
510 std::ifstream f(debug_filename_);
511 if (f.is_open())
512 std::cout << f.rdbuf();
513 }
514
515#ifdef WITH_VALGRIND
516 CALLGRIND_DUMP_STATS;
517#endif
518 }
519
520
521 /*
522 * ---------------------------------------------------------------------------
523 * Checkpointing, VTK output, and compute error:
524 * ---------------------------------------------------------------------------
525 */
526
527
528 template <typename Description, int dim, typename Number>
529 template <typename Callable>
531 StateVector &state_vector,
532 const std::string &base_name,
533 Number &t,
534 unsigned int &timer_cycle,
535 const Callable &prepare_compute_kernels)
536 {
537#ifdef DEBUG_OUTPUT
538 std::cout << "TimeLoop<dim, Number>::read_checkpoint()" << std::endl;
539#endif
540
541 AssertThrow(have_distributed_triangulation<dim>,
542 dealii::ExcMessage(
543 "read_checkpoint() is not implemented for "
544 "distributed::shared::Triangulation which we use in 1D"));
545
546 /*
547 * Initialize discretization, read in the mesh, and initialize everything:
548 */
549
550#if !DEAL_II_VERSION_GTE(9, 6, 0)
551 if constexpr (have_distributed_triangulation<dim>) {
552#endif
553 discretization_.refinement() = 0; /* do not refine */
554 discretization_.prepare(base_name);
555 discretization_.triangulation().load(base_name + "-checkpoint.mesh");
556#if !DEAL_II_VERSION_GTE(9, 6, 0)
557 }
558#endif
559
560 prepare_compute_kernels();
561
562 /*
563 * Read in and broadcast metadata:
564 */
565
566 std::string name = base_name + "-checkpoint";
567
568 unsigned int transfer_handle;
569 if (mpi_ensemble_.ensemble_rank() == 0) {
570 std::string meta = name + ".metadata";
571
572 std::ifstream file(meta, std::ios::binary);
573 boost::archive::binary_iarchive ia(file);
574 ia >> t >> timer_cycle >> transfer_handle;
575 }
576
577 int ierr;
578 if constexpr (std::is_same_v<Number, double>)
579 ierr = MPI_Bcast(
580 &t, 1, MPI_DOUBLE, 0, mpi_ensemble_.ensemble_communicator());
581 else
582 ierr =
583 MPI_Bcast(&t, 1, MPI_FLOAT, 0, mpi_ensemble_.ensemble_communicator());
584 AssertThrowMPI(ierr);
585
586 ierr = MPI_Bcast(&timer_cycle,
587 1,
588 MPI_UNSIGNED,
589 0,
590 mpi_ensemble_.ensemble_communicator());
591 AssertThrowMPI(ierr);
592
593 ierr = MPI_Bcast(&transfer_handle,
594 1,
595 MPI_UNSIGNED,
596 0,
597 mpi_ensemble_.ensemble_communicator());
598 AssertThrowMPI(ierr);
599
600 /* Now read in the state vector: */
601
602 Vectors::reinit_state_vector<Description>(state_vector, offline_data_);
603
604 solution_transfer_.set_handle(transfer_handle);
605 solution_transfer_.project(state_vector);
606 solution_transfer_.reset_handle();
607 }
608
609
610 template <typename Description, int dim, typename Number>
612 const StateVector &state_vector,
613 const std::string &base_name,
614 const Number &t,
615 const unsigned int &timer_cycle)
616 {
617#ifdef DEBUG_OUTPUT
618 std::cout << "TimeLoop<dim, Number>::write_checkpoint()" << std::endl;
619#endif
620
621 AssertThrow(have_distributed_triangulation<dim>,
622 dealii::ExcMessage(
623 "write_checkpoint() is not implemented for "
624 "distributed::shared::Triangulation which we use in 1D"));
625
626 /* We need hyperbolic_module.prepare_state_vector() prior to this call! */
627 solution_transfer_.prepare_projection(state_vector);
628 const auto transfer_handle = solution_transfer_.get_handle();
629 solution_transfer_.reset_handle();
630
631 std::string name = base_name + "-checkpoint";
632
633 if (mpi_ensemble_.ensemble_rank() == 0) {
634 for (const std::string suffix :
635 {".mesh", ".mesh_fixed.data", ".mesh.info", ".metadata"})
636 if (std::filesystem::exists(name + suffix))
637 std::filesystem::rename(name + suffix, name + suffix + "~");
638 }
639
640#if !DEAL_II_VERSION_GTE(9, 6, 0)
641 if constexpr (have_distributed_triangulation<dim>) {
642#endif
643 const auto &triangulation = discretization_.triangulation();
644 triangulation.save(name + ".mesh");
645#if !DEAL_II_VERSION_GTE(9, 6, 0)
646 }
647#endif
648
649 /*
650 * Now, write out metadata on rank 0:
651 */
652
653 if (mpi_ensemble_.ensemble_rank() == 0) {
654 std::string meta = name + ".metadata";
655 std::ofstream file(meta, std::ios::binary | std::ios::trunc);
656 boost::archive::binary_oarchive oa(file);
657 oa << t << timer_cycle << transfer_handle;
658 }
659
660 const int ierr = MPI_Barrier(mpi_ensemble_.ensemble_communicator());
661 AssertThrowMPI(ierr);
662 }
663
664
665 template <typename Description, int dim, typename Number>
666 template <typename Callable>
668 StateVector &state_vector, const Callable &prepare_compute_kernels)
669 {
670#ifdef DEBUG_OUTPUT
671 std::cout << "TimeLoop<dim, Number>::adapt_mesh_and_transfer_state_vector()"
672 << std::endl;
673#endif
674
675 AssertThrow(mpi_ensemble_.n_ensembles() == 1, dealii::ExcNotImplemented());
676
677 /*
678 * Mark cells for coarsening and refinement and set up triangulation:
679 */
680
681 auto &triangulation = discretization_.triangulation();
682 mesh_adaptor_.mark_cells_for_coarsening_and_refinement(triangulation);
683
684 triangulation.prepare_coarsening_and_refinement();
685
686 /* We need hyperbolic_module.prepare_state_vector() prior to this call! */
687 solution_transfer_.prepare_projection(state_vector);
688
689 /* Execute mesh adaptation and project old state to new state vector: */
690
691 triangulation.execute_coarsening_and_refinement();
692 prepare_compute_kernels();
693
694 Vectors::reinit_state_vector<Description>(state_vector, offline_data_);
695 solution_transfer_.project(state_vector);
696 solution_transfer_.reset_handle();
697 }
698
699
700 template <typename Description, int dim, typename Number>
701 void
703 const Number t)
704 {
705#ifdef DEBUG_OUTPUT
706 std::cout << "TimeLoop<dim, Number>::compute_error()" << std::endl;
707#endif
708
709 hyperbolic_module_.prepare_state_vector(state_vector, t);
710 if (!ParabolicSystem::is_identity)
711 parabolic_module_.prepare_state_vector(state_vector, t);
712
713 Vector<Number> difference_per_cell(
714 discretization_.triangulation().n_active_cells());
715
716 Number linf_norm = 0.;
717 Number l1_norm = 0;
718 Number l2_norm = 0;
719
720 const auto analytic_U =
721 initial_values_.get().interpolate_hyperbolic_vector(t);
722 const auto &U = std::get<0>(state_vector);
723
724 ScalarVector analytic_component;
725 ScalarVector error_component;
726 analytic_component.reinit(offline_data_.scalar_partitioner());
727 error_component.reinit(offline_data_.scalar_partitioner());
728
729 /* Loop over all selected components: */
730 for (const auto &entry : error_quantities_) {
731 const auto &names = View::component_names;
732 const auto pos = std::find(std::begin(names), std::end(names), entry);
733 if (pos == std::end(names)) {
734 AssertThrow(
735 false,
736 dealii::ExcMessage("Unknown component name »" + entry + "«"));
737 __builtin_trap();
738 }
739
740 const auto index = std::distance(std::begin(names), pos);
741
742 analytic_U.extract_component(analytic_component, index);
743
744 /* Compute norms of analytic solution: */
745
746 Number linf_norm_analytic = 0.;
747 Number l1_norm_analytic = 0.;
748 Number l2_norm_analytic = 0.;
749
750 if (error_normalize_) {
751 linf_norm_analytic =
752 Utilities::MPI::max(analytic_component.linfty_norm(),
753 mpi_ensemble_.ensemble_communicator());
754
755 VectorTools::integrate_difference(
756 offline_data_.dof_handler(),
757 analytic_component,
758 Functions::ZeroFunction<dim, Number>(),
759 difference_per_cell,
760 QGauss<dim>(3),
761 VectorTools::L1_norm);
762
763 l1_norm_analytic =
764 Utilities::MPI::sum(difference_per_cell.l1_norm(),
765 mpi_ensemble_.ensemble_communicator());
766
767 VectorTools::integrate_difference(
768 offline_data_.dof_handler(),
769 analytic_component,
770 Functions::ZeroFunction<dim, Number>(),
771 difference_per_cell,
772 QGauss<dim>(3),
773 VectorTools::L2_norm);
774
775 l2_norm_analytic = Number(std::sqrt(
776 Utilities::MPI::sum(std::pow(difference_per_cell.l2_norm(), 2),
777 mpi_ensemble_.ensemble_communicator())));
778 }
779
780 /* Compute norms of error: */
781
782 U.extract_component(error_component, index);
783 /* Populate constrained dofs due to periodicity: */
784 offline_data_.affine_constraints().distribute(error_component);
785 error_component.update_ghost_values();
786 error_component -= analytic_component;
787
788 const Number linf_norm_error = Utilities::MPI::max(
789 error_component.linfty_norm(), mpi_ensemble_.ensemble_communicator());
790
791 VectorTools::integrate_difference(offline_data_.dof_handler(),
792 error_component,
793 Functions::ZeroFunction<dim, Number>(),
794 difference_per_cell,
795 QGauss<dim>(3),
796 VectorTools::L1_norm);
797
798 const Number l1_norm_error = Utilities::MPI::sum(
799 difference_per_cell.l1_norm(), mpi_ensemble_.ensemble_communicator());
800
801 VectorTools::integrate_difference(offline_data_.dof_handler(),
802 error_component,
803 Functions::ZeroFunction<dim, Number>(),
804 difference_per_cell,
805 QGauss<dim>(3),
806 VectorTools::L2_norm);
807
808 const Number l2_norm_error = Number(std::sqrt(
809 Utilities::MPI::sum(std::pow(difference_per_cell.l2_norm(), 2),
810 mpi_ensemble_.ensemble_communicator())));
811
812 if (error_normalize_) {
813 linf_norm += linf_norm_error / linf_norm_analytic;
814 l1_norm += l1_norm_error / l1_norm_analytic;
815 l2_norm += l2_norm_error / l2_norm_analytic;
816 } else {
817 linf_norm += linf_norm_error;
818 l1_norm += l1_norm_error;
819 l2_norm += l2_norm_error;
820 }
821 }
822
823 if (mpi_ensemble_.ensemble_rank() != 0)
824 return;
825
826 /*
827 * Sum up over all participating MPI ranks. Note: we only perform this
828 * operation on "peer" ranks zero:
829 */
830
831 if (mpi_ensemble_.n_ensembles() > 1) {
832 linf_norm = Utilities::MPI::sum(
833 linf_norm, mpi_ensemble_.ensemble_leader_communicator());
834 l1_norm = Utilities::MPI::sum(
835 l1_norm, mpi_ensemble_.ensemble_leader_communicator());
836 l2_norm = Utilities::MPI::sum(
837 l2_norm, mpi_ensemble_.ensemble_leader_communicator());
838 }
839
840 if (mpi_ensemble_.world_rank() != 0)
841 return;
842
843 logfile_ << std::endl << "Computed errors:" << std::endl << std::endl;
844 logfile_ << std::setprecision(16);
845
846 std::string description =
847 error_normalize_ ? "Normalized consolidated" : "Consolidated";
848
849 logfile_ << description + " Linf, L1, and L2 errors at final time \n";
850 logfile_ << std::setprecision(16);
851 logfile_ << "#dofs = " << n_global_dofs_ << std::endl;
852 logfile_ << "t = " << t << std::endl;
853 logfile_ << "Linf = " << linf_norm << std::endl;
854 logfile_ << "L1 = " << l1_norm << std::endl;
855 logfile_ << "L2 = " << l2_norm << std::endl;
856
857 std::cout << description + " Linf, L1, and L2 errors at final time \n";
858 std::cout << std::setprecision(16);
859 std::cout << "#dofs = " << n_global_dofs_ << std::endl;
860 std::cout << "t = " << t << std::endl;
861 std::cout << "Linf = " << linf_norm << std::endl;
862 std::cout << "L1 = " << l1_norm << std::endl;
863 std::cout << "L2 = " << l2_norm << std::endl;
864 }
865
866
867 template <typename Description, int dim, typename Number>
869 const std::string &name,
870 const Number t,
871 const unsigned int cycle)
872 {
873#ifdef DEBUG_OUTPUT
874 std::cout << "TimeLoop<dim, Number>::output(t = " << t << ")" << std::endl;
875#endif
876
877 const bool do_full_output =
878 (cycle % timer_output_full_multiplier_ == 0) && enable_output_full_;
879 const bool do_levelsets =
880 (cycle % timer_output_levelsets_multiplier_ == 0) &&
881 enable_output_levelsets_;
882
883 /* There is nothing to do: */
884 if (!(do_full_output || do_levelsets))
885 return;
886
887 /* Prepare vectors for output: */
888
889 hyperbolic_module_.prepare_state_vector(state_vector, t);
890 if (!ParabolicSystem::is_identity)
891 parabolic_module_.prepare_state_vector(state_vector, t);
892
893 /* Data output: */
894
895 Scope scope(computing_timer_, "time step [X] - perform vtu output");
896 print_info("scheduling output");
897
898 postprocessor_.compute(state_vector);
899 /*
900 * Workaround: Manually reset bounds during the first output cycle
901 * (which is often just a uniform flow field) to obtain a better
902 * normailization:
903 */
904 if (cycle == 0)
905 postprocessor_.reset_bounds();
906
907 vtu_output_.schedule_output(
908 state_vector, name, t, cycle, do_full_output, do_levelsets);
909 }
910
911
912 /*
913 * ---------------------------------------------------------------------------
914 * Output and logging related functions:
915 * ---------------------------------------------------------------------------
916 */
917
918
919 template <typename Description, int dim, typename Number>
920 void
922 {
923 if (mpi_ensemble_.world_rank() != 0)
924 return;
925
926 /* Output commit and library information: */
927
929
930 /* Print run time parameters: */
931
932 stream << std::endl << "Run time parameters:" << std::endl << std::endl;
933 ParameterAcceptor::prm.print_parameters(
934 stream, ParameterHandler::OutputStyle::ShortPRM);
935 stream << std::endl;
936
937 /* Also print out parameters to a prm file: */
938
939 std::ofstream output(base_name_ + "-parameters.prm");
940 ParameterAcceptor::prm.print_parameters(output, ParameterHandler::ShortPRM);
941 }
942
943
944 template <typename Description, int dim, typename Number>
945 void
947 {
948 /*
949 * Fixme: this conversion to double is really not elegant. We should
950 * improve the Utilities::MPI::min_max_avg function in deal.II to
951 * handle different data types
952 */
953
954 // NOLINTBEGIN
955 std::vector<double> values = {
956 (double)offline_data_.n_export_indices(),
957 (double)offline_data_.n_locally_internal(),
958 (double)offline_data_.n_locally_owned(),
959 (double)offline_data_.n_locally_relevant(),
960 (double)offline_data_.n_export_indices() /
961 (double)offline_data_.n_locally_relevant(),
962 (double)offline_data_.n_locally_internal() /
963 (double)offline_data_.n_locally_relevant(),
964 (double)offline_data_.n_locally_owned() /
965 (double)offline_data_.n_locally_relevant()};
966 // NOLINTEND
967
968 const auto data =
969 Utilities::MPI::min_max_avg(values, mpi_ensemble_.world_communicator());
970
971 if (mpi_ensemble_.world_rank() != 0)
972 return;
973
974 std::ostringstream output;
975
976 unsigned int n =
977 dealii::Utilities::needed_digits(mpi_ensemble_.n_world_ranks());
978
979 const auto print_snippet = [&output, n](const std::string &name,
980 const auto &values) {
981 output << name << ": ";
982 // NOLINTBEGIN
983 output << std::setw(9) << (unsigned int)values.min //
984 << " [p" << std::setw(n) << values.min_index << "] " //
985 << std::setw(9) << (unsigned int)values.avg << " " //
986 << std::setw(9) << (unsigned int)values.max //
987 << " [p" << std::setw(n) << values.max_index << "]"; //
988 // NOLINTEND
989 };
990
991 const auto print_percentages = [&output, n](const auto &percentages) {
992 output << std::endl << " ";
993 output << " (" << std::setw(3) << std::setprecision(2)
994 << percentages.min * 100 << "% )"
995 << " [p" << std::setw(n) << percentages.min_index << "] "
996 << " (" << std::setw(3) << std::setprecision(2)
997 << percentages.avg * 100 << "% )"
998 << " "
999 << " (" << std::setw(3) << std::setprecision(2)
1000 << percentages.max * 100 << "% )"
1001 << " [p" << std::setw(n) << percentages.max_index << "]";
1002 };
1003
1004 output << std::endl << std::endl << "Partition: ";
1005 print_snippet("exp", data[0]);
1006 print_percentages(data[4]);
1007
1008 output << std::endl << " ";
1009 print_snippet("int", data[1]);
1010 print_percentages(data[5]);
1011
1012 output << std::endl << " ";
1013 print_snippet("own", data[2]);
1014 print_percentages(data[6]);
1015
1016 output << std::endl << " ";
1017 print_snippet("rel", data[3]);
1018
1019 stream << output.str() << std::endl;
1020 }
1021
1022
1023 template <typename Description, int dim, typename Number>
1025 {
1026 if (mpi_ensemble_.world_rank() != 0)
1027 return;
1028
1029 std::cout << "[INFO] " << header << std::endl;
1030 }
1031
1032
1033 template <typename Description, int dim, typename Number>
1034 void
1036 const std::string &secondary,
1037 std::ostream &stream)
1038 {
1039 if (mpi_ensemble_.world_rank() != 0)
1040 return;
1041
1042 const int header_size = header.size();
1043 const auto padded_header =
1044 std::string(std::max(0, 34 - header_size) / 2, ' ') + header +
1045 std::string(std::max(0, 35 - header_size) / 2, ' ');
1046
1047 const int secondary_size = secondary.size();
1048 const auto padded_secondary =
1049 std::string(std::max(0, 34 - secondary_size) / 2, ' ') + secondary +
1050 std::string(std::max(0, 35 - secondary_size) / 2, ' ');
1051
1052 /* clang-format off */
1053 stream << "\n";
1054 stream << " ####################################################\n";
1055 stream << " #########" << padded_header << "#########\n";
1056 stream << " #########" << padded_secondary << "#########\n";
1057 stream << " ####################################################\n";
1058 stream << std::endl;
1059 /* clang-format on */
1060 }
1061
1062
1063 template <typename Description, int dim, typename Number>
1065 unsigned int timer_cycle,
1066 Number last_checkpoint,
1067 std::ostream &stream,
1068 bool final_time)
1069 {
1070 static const std::string vectorization_name = [] {
1071 constexpr auto width = VectorizedArray<Number>::size();
1072
1073 std::string result;
1074 if (width == 1)
1075 result = "scalar ";
1076 else
1077 result = std::to_string(width * 8 * sizeof(Number)) + " bit packed ";
1078
1079 if constexpr (std::is_same_v<Number, double>)
1080 return result + "double";
1081 else if constexpr (std::is_same_v<Number, float>)
1082 return result + "float";
1083 else
1084 __builtin_trap();
1085 }();
1086
1087 stream << "Information: (HYP) " << hyperbolic_system_.get().problem_name;
1088 if constexpr (!ParabolicSystem::is_identity) {
1089 stream << "\n (PAR) " << parabolic_system_.get().problem_name;
1090 }
1091 stream << "\n [" << base_name_ << "] ";
1092 if (mpi_ensemble_.n_ensembles() > 1) {
1093 stream << mpi_ensemble_.n_ensembles() << " ensembles ";
1094 }
1095 stream << "with " //
1096 << n_global_dofs_ << " Qdofs on " //
1097 << mpi_ensemble_.n_world_ranks() << " ranks / " //
1098#ifdef WITH_OPENMP
1099 << MultithreadInfo::n_threads() << " threads <" //
1100#else
1101 << "[openmp disabled] <" //
1102#endif
1103 << vectorization_name << ">\n";
1104
1105 stream << " Last output cycle " //
1106 << timer_cycle - 1 //
1107 << " at t = " << timer_granularity_ * (timer_cycle - 1) //
1108 << " [ log ";
1109
1110 if (enable_output_full_)
1111 stream << "full ";
1112 if (enable_output_levelsets_)
1113 stream << "levelsets ";
1114 if (enable_compute_quantities_)
1115 stream << "quantities ";
1116
1117 stream << "]\n";
1118
1119 if (checkpoint_update_interval_ != Number(0.)) {
1120 const auto wall_time =
1121 Utilities::MPI::min_max_avg(computing_timer_["time loop"].wall_time(),
1122 mpi_ensemble_.world_communicator());
1123
1124 if (final_time) {
1125 stream << " Last checkpoint at FINAL TIME\n";
1126 } else {
1127 stream << " Last checkpoint at wall time " //
1128 << std::setprecision(2) << std::fixed << last_checkpoint //
1129 << "s (" << std::setprecision(0)
1130 << std::max(0., wall_time.max - last_checkpoint)
1131 << "s ago, interval " << checkpoint_update_interval_ << "s)\n";
1132 }
1133 }
1134 }
1135
1136
1137 template <typename Description, int dim, typename Number>
1139 std::ostream &stream)
1140 {
1141 Utilities::System::MemoryStats stats;
1142 Utilities::System::get_memory_stats(stats);
1143
1144 Utilities::MPI::MinMaxAvg data = Utilities::MPI::min_max_avg(
1145 stats.VmRSS / 1024., mpi_ensemble_.world_communicator());
1146
1147 if (mpi_ensemble_.world_rank() != 0)
1148 return;
1149
1150 std::ostringstream output;
1151
1152 unsigned int n =
1153 dealii::Utilities::needed_digits(mpi_ensemble_.n_world_ranks());
1154
1155 output << "\nMemory: [MiB]" //
1156 << std::setw(8) << data.min //
1157 << " [p" << std::setw(n) << data.min_index << "] " //
1158 << std::setw(8) << data.avg << " " //
1159 << std::setw(8) << data.max //
1160 << " [p" << std::setw(n) << data.max_index << "]"; //
1161
1162 stream << output.str() << std::endl;
1163 }
1164
1165
1166 template <typename Description, int dim, typename Number>
1168 {
1169 std::vector<std::ostringstream> output(computing_timer_.size());
1170
1171 const auto equalize = [&]() {
1172 const auto ptr =
1173 std::max_element(output.begin(),
1174 output.end(),
1175 [](const auto &left, const auto &right) {
1176 return left.str().length() < right.str().length();
1177 });
1178 const auto length = ptr->str().length();
1179 for (auto &it : output)
1180 it << std::string(length - it.str().length() + 1, ' ');
1181 };
1182
1183 const auto print_wall_time = [&](auto &timer, auto &stream) {
1184 const auto wall_time = Utilities::MPI::min_max_avg(
1185 timer.wall_time(), mpi_ensemble_.world_communicator());
1186
1187 constexpr auto eps = std::numeric_limits<double>::epsilon();
1188 /*
1189 * Cut off at 99.9% to avoid silly percentages cluttering up the
1190 * output.
1191 */
1192 const auto skew_negative = std::max(
1193 100. * (wall_time.min - wall_time.avg) / wall_time.avg - eps, -99.9);
1194 const auto skew_positive = std::min(
1195 100. * (wall_time.max - wall_time.avg) / wall_time.avg + eps, 99.9);
1196
1197 stream << std::setprecision(2) << std::fixed << std::setw(9)
1198 << wall_time.avg << "s [sk: " << std::setprecision(1)
1199 << std::setw(5) << std::fixed << skew_negative << "%/"
1200 << std::setw(4) << std::fixed << skew_positive << "%]";
1201 unsigned int n =
1202 dealii::Utilities::needed_digits(mpi_ensemble_.n_world_ranks());
1203 stream << " [p" << std::setw(n) << wall_time.min_index << "/"
1204 << wall_time.max_index << "]";
1205 };
1206
1207 const auto cpu_time_statistics =
1208 Utilities::MPI::min_max_avg(computing_timer_["time loop"].cpu_time(),
1209 mpi_ensemble_.world_communicator());
1210 const double total_cpu_time = cpu_time_statistics.sum;
1211
1212 const auto print_cpu_time =
1213 [&](auto &timer, auto &stream, bool percentage) {
1214 const auto cpu_time = Utilities::MPI::min_max_avg(
1215 timer.cpu_time(), mpi_ensemble_.world_communicator());
1216
1217 stream << std::setprecision(2) << std::fixed << std::setw(12)
1218 << cpu_time.sum << "s ";
1219
1220 if (percentage)
1221 stream << "(" << std::setprecision(1) << std::setw(4)
1222 << 100. * cpu_time.sum / total_cpu_time << "%)";
1223 };
1224
1225 auto jt = output.begin();
1226 for (auto &it : computing_timer_)
1227 *jt++ << " " << it.first;
1228 equalize();
1229
1230 jt = output.begin();
1231 for (auto &it : computing_timer_)
1232 print_wall_time(it.second, *jt++);
1233 equalize();
1234
1235 jt = output.begin();
1236 bool compute_percentages = false;
1237 for (auto &it : computing_timer_) {
1238 print_cpu_time(it.second, *jt++, compute_percentages);
1239 if (it.first.starts_with("time loop"))
1240 compute_percentages = true;
1241 }
1242 equalize();
1243
1244 if (mpi_ensemble_.world_rank() != 0)
1245 return;
1246
1247 stream << std::endl << "Timer statistics:\n";
1248 for (auto &it : output)
1249 stream << it.str() << std::endl;
1250 }
1251
1252
1253 template <typename Description, int dim, typename Number>
1255 unsigned int cycle, Number t, std::ostream &stream, bool final_time)
1256 {
1257 /*
1258 * Fixme: The global state kept in this function should be refactored
1259 * into its own class object.
1260 */
1261 static struct Data {
1262 unsigned int cycle = 0;
1263 double t = 0.;
1264 double cpu_time_sum = 0.;
1265 double cpu_time_avg = 0.;
1266 double cpu_time_min = 0.;
1267 double cpu_time_max = 0.;
1268 double wall_time = 0.;
1269 } previous, current;
1270
1271 static double time_per_second_exp = 0.;
1272
1273 /* Update statistics: */
1274
1275 {
1276 previous = current;
1277
1278 current.cycle = cycle;
1279 current.t = t;
1280
1281 const auto wall_time_statistics =
1282 Utilities::MPI::min_max_avg(computing_timer_["time loop"].wall_time(),
1283 mpi_ensemble_.world_communicator());
1284 current.wall_time = wall_time_statistics.max;
1285
1286 const auto cpu_time_statistics =
1287 Utilities::MPI::min_max_avg(computing_timer_["time loop"].cpu_time(),
1288 mpi_ensemble_.world_communicator());
1289 current.cpu_time_sum = cpu_time_statistics.sum;
1290 current.cpu_time_avg = cpu_time_statistics.avg;
1291 current.cpu_time_min = cpu_time_statistics.min;
1292 current.cpu_time_max = cpu_time_statistics.max;
1293 }
1294
1295 if (final_time)
1296 previous = Data();
1297
1298 /* Take averages: */
1299
1300 double delta_cycles = current.cycle - previous.cycle;
1301 const double cycles_per_second =
1302 delta_cycles / (current.wall_time - previous.wall_time);
1303
1304 const auto efficiency = time_integrator_.efficiency();
1305 const auto n_dofs = static_cast<double>(n_global_dofs_);
1306
1307 const double wall_m_dofs_per_sec =
1308 delta_cycles * n_dofs / 1.e6 /
1309 (current.wall_time - previous.wall_time) * efficiency;
1310
1311 double cpu_m_dofs_per_sec = delta_cycles * n_dofs / 1.e6 /
1312 (current.cpu_time_sum - previous.cpu_time_sum) *
1313 efficiency;
1314#ifdef WITH_OPENMP
1315 if (terminal_show_rank_throughput_)
1316 cpu_m_dofs_per_sec *= MultithreadInfo::n_threads();
1317#endif
1318
1319 double cpu_time_skew = (current.cpu_time_max - current.cpu_time_min - //
1320 previous.cpu_time_max + previous.cpu_time_min) /
1321 delta_cycles;
1322 /* avoid printing small negative numbers: */
1323 cpu_time_skew = std::max(0., cpu_time_skew);
1324
1325 const double cpu_time_skew_percentage =
1326 cpu_time_skew * delta_cycles /
1327 (current.cpu_time_avg - previous.cpu_time_avg);
1328
1329 const double delta_time =
1330 (current.t - previous.t) / (current.cycle - previous.cycle);
1331 const double time_per_second =
1332 (current.t - previous.t) / (current.wall_time - previous.wall_time);
1333
1334 /* Print Jean-Luc and Martin metrics: */
1335
1336 std::ostringstream output;
1337
1338 /* clang-format off */
1339 output << std::endl;
1340
1341 output << "Throughput:\n "
1342 << (terminal_show_rank_throughput_? "RANK: " : "CPU : ")
1343 << std::setprecision(4) << std::fixed << cpu_m_dofs_per_sec
1344 << " MQ/s ("
1345 << std::scientific << 1. / cpu_m_dofs_per_sec * 1.e-6
1346 << " s/Qdof/substep)" << std::endl;
1347
1348 output << " [cpu time skew: "
1349 << std::setprecision(2) << std::scientific << cpu_time_skew
1350 << "s/cycle ("
1351 << std::setprecision(1) << std::setw(4) << std::setfill(' ') << std::fixed
1352 << 100. * cpu_time_skew_percentage
1353 << "%)]" << std::endl;
1354
1355 output << " WALL: "
1356 << std::setprecision(4) << std::fixed << wall_m_dofs_per_sec
1357 << " MQ/s ("
1358 << std::scientific << 1. / wall_m_dofs_per_sec * 1.e-6
1359 << " s/Qdof/substep) ("
1360 << std::setprecision(2) << std::fixed << cycles_per_second
1361 << " cycles/s)" << std::endl;
1362
1363 const auto &scheme = time_integrator_.time_stepping_scheme();
1364 output << " [ "
1365 << Patterns::Tools::Convert<TimeSteppingScheme>::to_string(scheme)
1366 << " with CFL = "
1367 << std::setprecision(2) << std::fixed << hyperbolic_module_.cfl()
1368 << " ("
1369 << std::setprecision(0) << std::fixed << hyperbolic_module_.n_restarts()
1370 << "/"
1371 << std::setprecision(0) << std::fixed << parabolic_module_.n_restarts()
1372 << " rsts) ("
1373 << std::setprecision(0) << std::fixed << hyperbolic_module_.n_warnings()
1374 << "/"
1375 << std::setprecision(0) << std::fixed << parabolic_module_.n_warnings()
1376 << " warn) ("
1377 << std::setprecision(0) << std::fixed << hyperbolic_module_.n_corrections()
1378 << "/"
1379 << std::setprecision(0) << std::fixed << parabolic_module_.n_corrections()
1380 << " corr) ]" << std::endl;
1381
1382 if constexpr (!ParabolicSystem::is_identity)
1383 parabolic_module_.print_solver_statistics(output);
1384
1385 output << " [ dt = "
1386 << std::scientific << std::setprecision(2) << delta_time
1387 << " ( "
1388 << time_per_second
1389 << " dt/s) ]" << std::endl;
1390 /* clang-format on */
1391
1392 /* And print an ETA: */
1393
1394 time_per_second_exp = 0.8 * time_per_second_exp + 0.2 * time_per_second;
1395 auto eta = static_cast<unsigned int>(std::max(t_final_ - t, Number(0.)) /
1396 time_per_second_exp);
1397
1398 output << "\n ETA : ";
1399
1400 const unsigned int days = eta / (24 * 3600);
1401 if (days > 0) {
1402 output << days << " d ";
1403 eta %= 24 * 3600;
1404 }
1405
1406 const unsigned int hours = eta / 3600;
1407 if (hours > 0) {
1408 output << hours << " h ";
1409 eta %= 3600;
1410 }
1411
1412 const unsigned int minutes = eta / 60;
1413 output << minutes << " min";
1414
1415 output << " (terminal update every " //
1416 << std::setprecision(2) << std::fixed << terminal_update_interval_
1417 << "s)";
1418
1419 if (mpi_ensemble_.world_rank() != 0)
1420 return;
1421
1422 stream << output.str() << std::endl;
1423 }
1424
1425
1426 template <typename Description, int dim, typename Number>
1428 unsigned int cycle,
1429 Number t,
1430 unsigned int timer_cycle,
1431 Number last_checkpoint,
1432 bool write_to_logfile,
1433 bool final_time)
1434 {
1435 std::ostringstream output;
1436
1437 /* Print header: */
1438
1439 std::ostringstream primary;
1440 if (final_time) {
1441 primary << "FINAL (cycle " << Utilities::int_to_string(cycle, 6) << ")";
1442 } else {
1443 primary << "Cycle " << Utilities::int_to_string(cycle, 6) //
1444 << " (" << std::fixed << std::setprecision(1) //
1445 << t / t_final_ * 100 << "%)";
1446 }
1447
1448 std::ostringstream secondary;
1449 secondary << "at time t = " << std::setprecision(8) << std::fixed << t;
1450
1451 print_head(primary.str(), secondary.str(), output);
1452
1453 /* Print information and statistics: */
1454
1455 print_information(timer_cycle, last_checkpoint, output, final_time);
1456 print_memory_statistics(output);
1457 print_timers(output);
1458 print_throughput(cycle, t, output, final_time);
1459
1460 /* Only output on rank 0: */
1461 if (mpi_ensemble_.world_rank() != 0)
1462 return;
1463
1464#ifndef DEBUG_OUTPUT
1465 std::cout << "\033[2J\033[H";
1466#endif
1467 std::cout << output.str() << std::flush;
1468
1469 if (write_to_logfile) {
1470 logfile_ << "\n" << output.str() << std::flush;
1471 }
1472 }
1473} // namespace ryujin
void write_checkpoint(const StateVector &state_vector, const std::string &base_name, const Number &t, const unsigned int &output_cycle)
Vectors::ScalarVector< Number > ScalarVector
Definition: time_loop.h:63
void print_timers(std::ostream &stream)
void output(StateVector &state_vector, const std::string &name, const Number t, const unsigned int cycle)
void print_memory_statistics(std::ostream &stream)
void print_mpi_partition(std::ostream &stream)
void print_parameters(std::ostream &stream)
void compute_error(StateVector &state_vector, Number t)
void read_checkpoint(StateVector &state_vector, const std::string &base_name, Number &t, unsigned int &output_cycle, const Callable &prepare_compute_kernels)
TimeLoop(const MPI_Comm &mpi_comm)
void print_cycle_statistics(unsigned int cycle, Number t, unsigned int output_cycle, Number last_checkpoint, bool write_to_logfile=false, bool final_time=false)
typename View::StateVector StateVector
Definition: time_loop.h:61
void print_head(const std::string &header, const std::string &secondary, std::ostream &stream)
void print_throughput(unsigned int cycle, Number t, std::ostream &stream, bool final_time=false)
void print_information(unsigned int output_cycle, Number last_checkpoint, std::ostream &stream, bool final_time=false)
void print_info(const std::string &header)
void adapt_mesh_and_transfer_state_vector(StateVector &state_vector, const Callable &prepare_compute_kernels)
T pow(const T x, const T b)
void print_revision_and_version(std::ostream &stream)