13#include <deal.II/base/logstream.h>
14#include <deal.II/base/work_stream.h>
15#include <deal.II/numerics/vector_tools.h>
16#include <deal.II/numerics/vector_tools.templates.h>
22using namespace dealii;
26 template <
typename Description,
int dim,
typename Number>
28 : ParameterAcceptor(
"/A - TimeLoop")
29 , mpi_ensemble_(mpi_comm)
30 , hyperbolic_system_(mpi_ensemble_,
"/B - Equation")
31 , parabolic_system_(mpi_ensemble_,
"/B - Equation")
32 , discretization_(mpi_ensemble_,
"/C - Discretization")
33 , offline_data_(mpi_ensemble_, discretization_,
"/D - OfflineData")
34 , initial_values_(mpi_ensemble_,
40 , hyperbolic_module_(mpi_ensemble_,
45 "/F - HyperbolicModule")
46 , parabolic_module_(mpi_ensemble_,
52 "/G - ParabolicModule")
53 , time_integrator_(mpi_ensemble_,
57 "/H - TimeIntegrator")
58 , mesh_adaptor_(mpi_ensemble_,
62 hyperbolic_module_.initial_precomputed(),
63 hyperbolic_module_.alpha(),
66 mpi_ensemble_, offline_data_, hyperbolic_system_, parabolic_system_)
67 , postprocessor_(mpi_ensemble_,
72 , vtu_output_(mpi_ensemble_,
77 hyperbolic_module_.initial_precomputed(),
78 hyperbolic_module_.alpha(),
80 , quantities_(mpi_ensemble_,
87 add_parameter(
"basename", base_name_,
"Base name for all output files");
89 t_final_ = Number(5.);
90 add_parameter(
"final time", t_final_,
"Final time");
92 enforce_t_final_ =
false;
93 add_parameter(
"enforce final time",
95 "Boolean indicating whether the final time should be "
96 "enforced strictly. If set to true the last time step is "
97 "shortened so that the simulation ends precisely at t_final");
99 timer_granularity_ = Number(0.01);
100 add_parameter(
"timer granularity",
102 "The timer granularity specifies the time interval after "
103 "which compute, output, postprocessing, and mesh adaptation "
104 "routines are run. This \"baseline tick\" is further "
105 "modified by the corresponding \"*_multiplier\" options");
107 enable_checkpointing_ =
false;
109 "enable checkpointing",
110 enable_checkpointing_,
111 "Write out checkpoints to resume an interrupted computation at timer "
112 "granularity intervals. The frequency is determined by \"timer "
113 "granularity\" and \"timer checkpoint multiplier\"");
115 enable_output_full_ =
false;
116 add_parameter(
"enable output full",
118 "Write out full pvtu records. The frequency is determined by "
119 "\"timer granularity\" and \"timer output full multiplier\"");
121 enable_output_levelsets_ =
false;
123 "enable output levelsets",
124 enable_output_levelsets_,
125 "Write out levelsets pvtu records. The frequency is determined by "
126 "\"timer granularity\" and \"timer output levelsets multiplier\"");
128 enable_compute_error_ =
false;
129 add_parameter(
"enable compute error",
130 enable_compute_error_,
131 "Flag to control whether we compute the Linfty Linf_norm of "
132 "the difference to an analytic solution. Implemented only "
133 "for certain initial state configurations.");
135 enable_compute_quantities_ =
false;
137 "enable compute quantities",
138 enable_compute_quantities_,
139 "Flag to control whether we compute quantities of interest. The "
140 "frequency how often quantities are logged is determined by \"timer "
141 "granularity\" and \"timer compute quantities multiplier\"");
143 enable_mesh_adaptivity_ =
false;
145 "enable mesh adaptivity",
146 enable_mesh_adaptivity_,
147 "Flag to control whether we use an adaptive mesh refinement strategy. "
148 "The frequency how often we query MeshAdaptor::analyze() for deciding "
149 "on adapting the mesh is determined by \"timer granularity\" and "
150 "\"timer mesh refinement multiplier\"");
152 timer_checkpoint_multiplier_ = 1;
153 add_parameter(
"timer checkpoint multiplier",
154 timer_checkpoint_multiplier_,
155 "Multiplicative modifier applied to \"timer granularity\" "
156 "that determines the checkpointing granularity");
158 timer_output_full_multiplier_ = 1;
159 add_parameter(
"timer output full multiplier",
160 timer_output_full_multiplier_,
161 "Multiplicative modifier applied to \"timer granularity\" "
162 "that determines the full pvtu writeout granularity");
164 timer_output_levelsets_multiplier_ = 1;
165 add_parameter(
"timer output levelsets multiplier",
166 timer_output_levelsets_multiplier_,
167 "Multiplicative modifier applied to \"timer granularity\" "
168 "that determines the levelsets pvtu writeout granularity");
170 timer_compute_quantities_multiplier_ = 1;
172 "timer compute quantities multiplier",
173 timer_compute_quantities_multiplier_,
174 "Multiplicative modifier applied to \"timer granularity\" that "
175 "determines the writeout granularity for quantities of interest");
177 std::copy(std::begin(View::component_names),
178 std::end(View::component_names),
179 std::back_inserter(error_quantities_));
181 add_parameter(
"error quantities",
183 "List of conserved quantities used in the computation of the "
186 error_normalize_ =
true;
187 add_parameter(
"error normalize",
189 "Flag to control whether the error should be normalized by "
190 "the corresponding norm of the analytic solution.");
193 add_parameter(
"resume", resume_,
"Resume an interrupted computation");
195 resume_at_time_zero_ =
false;
196 add_parameter(
"resume at time zero",
197 resume_at_time_zero_,
198 "Resume from the latest checkpoint but set the time to t=0.");
200 terminal_update_interval_ = 5;
201 add_parameter(
"terminal update interval",
202 terminal_update_interval_,
203 "Number of seconds after which output statistics are "
204 "recomputed and printed on the terminal");
206 terminal_show_rank_throughput_ =
true;
207 add_parameter(
"terminal show rank throughput",
208 terminal_show_rank_throughput_,
209 "If set to true an average per rank throughput is computed "
210 "by dividing the total consumed CPU time (per rank) by the "
211 "number of threads (per rank). If set to false then a plain "
212 "average per thread \"CPU\" throughput value is computed by "
213 "using the umodified total accumulated CPU time.");
215 debug_filename_ =
"";
216 add_parameter(
"debug filename",
218 "If set to a nonempty string then we output the contents of "
219 "this file at the end. This is mainly useful in the "
220 "testsuite to output files we wish to compare");
224 template <
typename Description,
int dim,
typename Number>
228 std::cout <<
"TimeLoop<dim, Number>::run()" << std::endl;
232 base_name_ensemble_ = base_name_;
233 if (mpi_ensemble_.n_ensembles() > 1) {
234 print_info(
"setting up MPI ensemble");
235 base_name_ensemble_ +=
"-ensemble_" + dealii::Utilities::int_to_string(
236 mpi_ensemble_.ensemble(),
237 mpi_ensemble_.n_ensembles());
243 if (mpi_ensemble_.world_rank() == 0)
244 logfile_.open(base_name_ +
".log");
246 print_parameters(logfile_);
253 unsigned int timer_cycle = 0;
257 const auto prepare_compute_kernels = [&]() {
258 print_info(
"preparing compute kernels");
260 unsigned int n_parabolic_state_vectors =
261 parabolic_system_.get().n_parabolic_state_vectors();
263 offline_data_.prepare(
264 problem_dimension, n_precomputed_values, n_parabolic_state_vectors);
266 hyperbolic_module_.prepare();
267 parabolic_module_.prepare();
268 time_integrator_.prepare();
269 mesh_adaptor_.prepare( t);
270 postprocessor_.prepare();
271 vtu_output_.prepare();
272 quantities_.prepare(base_name_ensemble_);
273 print_mpi_partition(logfile_);
275 if (mpi_ensemble_.ensemble_rank() == 0)
276 n_global_dofs_ = dealii::Utilities::MPI::sum(
277 offline_data_.dof_handler().n_dofs(),
278 mpi_ensemble_.ensemble_leader_communicator());
282 Scope scope(computing_timer_,
"(re)initialize data structures");
283 print_info(
"initializing data structures");
286 print_info(
"resume: reading mesh and loading state vector");
288 read_checkpoint(state_vector,
292 prepare_compute_kernels);
294 if (resume_at_time_zero_) {
301 print_info(
"creating mesh and interpolating initial values");
303 discretization_.prepare(base_name_ensemble_);
305 prepare_compute_kernels();
307 Vectors::reinit_state_vector<Description>(state_vector, offline_data_);
308 std::get<0>(state_vector) =
309 initial_values_.get().interpolate_hyperbolic_vector();
317 Vectors::debug_poison_constrained_dofs<Description>(state_vector,
319 Vectors::debug_poison_precomputed_values<Description>(state_vector,
322 unsigned int cycle = 1;
323 Number last_terminal_output = (terminal_update_interval_ == Number(0.)
324 ? std::numeric_limits<Number>::max()
325 : std::numeric_limits<Number>::lowest());
331 print_info(
"entering main loop");
332 computing_timer_[
"time loop"].start();
334 constexpr Number relax =
335 Number(1.) - Number(10.) * std::numeric_limits<Number>::epsilon();
340 std::cout <<
"\n\n### cycle = " << cycle <<
" ###\n\n" << std::endl;
345 if (enable_compute_quantities_) {
346 Scope scope(computing_timer_,
347 "time step [X] - accumulate quantities");
348 quantities_.accumulate(state_vector, t);
353 if (t >= relax * timer_cycle * timer_granularity_) {
354 if (enable_compute_error_) {
362 Scope scope(computing_timer_,
363 "time step [X] - interpolate analytic solution");
364 Vectors::reinit_state_vector<Description>(analytic, offline_data_);
365 std::get<0>(analytic) =
366 initial_values_.get().interpolate_hyperbolic_vector(t);
375 base_name_ensemble_ +
"-analytic_solution",
380 output(state_vector, base_name_ensemble_ +
"-solution", t, timer_cycle);
382 if (enable_compute_quantities_ &&
383 (timer_cycle % timer_compute_quantities_multiplier_ == 0)) {
384 Scope scope(computing_timer_,
385 "time step [X] - write out quantities");
386 quantities_.write_out(state_vector, t, timer_cycle);
394 if (t >= relax * t_final_)
399 if (enable_mesh_adaptivity_) {
401 Scope scope(computing_timer_,
402 "time step [X] - analyze for mesh adaptation");
404 mesh_adaptor_.analyze(state_vector, t, cycle);
407 if (mesh_adaptor_.need_mesh_adaptation()) {
408 Scope scope(computing_timer_,
"(re)initialize data structures");
409 print_info(
"performing mesh adaptation");
411 hyperbolic_module_.prepare_state_vector(state_vector, t);
412 if (!ParabolicSystem::is_identity)
413 parabolic_module_.prepare_state_vector(state_vector, t);
414 adapt_mesh_and_transfer_state_vector(state_vector,
415 prepare_compute_kernels);
421 const auto tau = time_integrator_.step(
425 ? std::min(t_final_, timer_cycle * timer_granularity_)
426 : std::numeric_limits<Number>::max());
431 if (terminal_update_interval_ != Number(0.)) {
434 const bool write_to_log_file =
435 (t >= relax * timer_cycle * timer_granularity_);
438 const auto wall_time = computing_timer_[
"time loop"].wall_time();
439 int update_terminal =
440 (wall_time >= last_terminal_output + terminal_update_interval_);
443 const auto ierr = MPI_Bcast(&update_terminal,
447 mpi_ensemble_.world_communicator());
448 AssertThrowMPI(ierr);
450 if (write_to_log_file || update_terminal) {
451 print_cycle_statistics(
452 cycle, t, timer_cycle, write_to_log_file);
453 last_terminal_output = wall_time;
461 computing_timer_[
"time loop"].stop();
463 if (terminal_update_interval_ != Number(0.)) {
465 print_cycle_statistics(
466 cycle, t, timer_cycle,
true,
true);
469 if (enable_compute_error_) {
471 compute_error(state_vector, t);
474 if (mpi_ensemble_.world_rank() == 0 && debug_filename_ !=
"") {
475 std::ifstream f(debug_filename_);
477 std::cout << f.rdbuf();
481 CALLGRIND_DUMP_STATS;
486 template <
typename Description,
int dim,
typename Number>
487 template <
typename Callable>
490 const std::string &base_name,
492 unsigned int &output_cycle,
493 const Callable &prepare_compute_kernels)
496 std::cout <<
"TimeLoop<dim, Number>::read_checkpoint()" << std::endl;
499 AssertThrow(have_distributed_triangulation<dim>,
501 "read_checkpoint() is not implemented for "
502 "distributed::shared::Triangulation which we use in 1D"));
508#if !DEAL_II_VERSION_GTE(9, 6, 0)
509 if constexpr (have_distributed_triangulation<dim>) {
511 discretization_.refinement() = 0;
512 discretization_.prepare(base_name);
513 discretization_.triangulation().load(base_name +
"-checkpoint.mesh");
514#if !DEAL_II_VERSION_GTE(9, 6, 0)
518 prepare_compute_kernels();
524 std::string name = base_name +
"-checkpoint";
526 unsigned int transfer_handle;
527 if (mpi_ensemble_.ensemble_rank() == 0) {
528 std::string meta = name +
".metadata";
530 std::ifstream file(meta, std::ios::binary);
531 boost::archive::binary_iarchive ia(file);
532 ia >> t >> output_cycle >> transfer_handle;
536 if constexpr (std::is_same_v<Number, double>)
538 &t, 1, MPI_DOUBLE, 0, mpi_ensemble_.ensemble_communicator());
541 MPI_Bcast(&t, 1, MPI_FLOAT, 0, mpi_ensemble_.ensemble_communicator());
542 AssertThrowMPI(ierr);
544 ierr = MPI_Bcast(&output_cycle,
548 mpi_ensemble_.ensemble_communicator());
549 AssertThrowMPI(ierr);
551 ierr = MPI_Bcast(&transfer_handle,
555 mpi_ensemble_.ensemble_communicator());
556 AssertThrowMPI(ierr);
560 Vectors::reinit_state_vector<Description>(state_vector, offline_data_);
562 solution_transfer_.set_handle(transfer_handle);
563 solution_transfer_.project(state_vector);
564 solution_transfer_.reset_handle();
568 template <
typename Description,
int dim,
typename Number>
571 const std::string &base_name,
573 const unsigned int &output_cycle)
576 std::cout <<
"TimeLoop<dim, Number>::write_checkpoint()" << std::endl;
579 AssertThrow(have_distributed_triangulation<dim>,
581 "write_checkpoint() is not implemented for "
582 "distributed::shared::Triangulation which we use in 1D"));
585 solution_transfer_.prepare_projection(state_vector);
586 const auto transfer_handle = solution_transfer_.get_handle();
587 solution_transfer_.reset_handle();
589 std::string name = base_name +
"-checkpoint";
591 if (mpi_ensemble_.ensemble_rank() == 0) {
592 for (
const std::string suffix :
593 {
".mesh",
".mesh_fixed.data",
".mesh.info",
".metadata"})
594 if (std::filesystem::exists(name + suffix))
595 std::filesystem::rename(name + suffix, name + suffix +
"~");
598#if !DEAL_II_VERSION_GTE(9, 6, 0)
599 if constexpr (have_distributed_triangulation<dim>) {
601 const auto &triangulation = discretization_.triangulation();
602 triangulation.save(name +
".mesh");
603#if !DEAL_II_VERSION_GTE(9, 6, 0)
611 if (mpi_ensemble_.ensemble_rank() == 0) {
612 std::string meta = name +
".metadata";
613 std::ofstream file(meta, std::ios::binary | std::ios::trunc);
614 boost::archive::binary_oarchive oa(file);
615 oa << t << output_cycle << transfer_handle;
618 const int ierr = MPI_Barrier(mpi_ensemble_.ensemble_communicator());
619 AssertThrowMPI(ierr);
623 template <
typename Description,
int dim,
typename Number>
624 template <
typename Callable>
626 StateVector &state_vector,
const Callable &prepare_compute_kernels)
629 std::cout <<
"TimeLoop<dim, Number>::adapt_mesh_and_transfer_state_vector()"
633 AssertThrow(mpi_ensemble_.n_ensembles() == 1, dealii::ExcNotImplemented());
639 auto &triangulation = discretization_.triangulation();
640 mesh_adaptor_.mark_cells_for_coarsening_and_refinement(triangulation);
642 triangulation.prepare_coarsening_and_refinement();
645 solution_transfer_.prepare_projection(state_vector);
649 triangulation.execute_coarsening_and_refinement();
650 prepare_compute_kernels();
652 Vectors::reinit_state_vector<Description>(state_vector, offline_data_);
653 solution_transfer_.project(state_vector);
654 solution_transfer_.reset_handle();
658 template <
typename Description,
int dim,
typename Number>
664 std::cout <<
"TimeLoop<dim, Number>::compute_error()" << std::endl;
667 hyperbolic_module_.prepare_state_vector(state_vector, t);
668 if (!ParabolicSystem::is_identity)
669 parabolic_module_.prepare_state_vector(state_vector, t);
671 Vector<Number> difference_per_cell(
672 discretization_.triangulation().n_active_cells());
674 Number linf_norm = 0.;
678 const auto analytic_U =
679 initial_values_.get().interpolate_hyperbolic_vector(t);
680 const auto &U = std::get<0>(state_vector);
684 analytic_component.reinit(offline_data_.scalar_partitioner());
685 error_component.reinit(offline_data_.scalar_partitioner());
688 for (
const auto &entry : error_quantities_) {
689 const auto &names = View::component_names;
690 const auto pos = std::find(std::begin(names), std::end(names), entry);
691 if (pos == std::end(names)) {
694 dealii::ExcMessage(
"Unknown component name »" + entry +
"«"));
698 const auto index = std::distance(std::begin(names), pos);
700 analytic_U.extract_component(analytic_component, index);
704 Number linf_norm_analytic = 0.;
705 Number l1_norm_analytic = 0.;
706 Number l2_norm_analytic = 0.;
708 if (error_normalize_) {
710 Utilities::MPI::max(analytic_component.linfty_norm(),
711 mpi_ensemble_.ensemble_communicator());
713 VectorTools::integrate_difference(
714 offline_data_.dof_handler(),
716 Functions::ZeroFunction<dim, Number>(),
719 VectorTools::L1_norm);
722 Utilities::MPI::sum(difference_per_cell.l1_norm(),
723 mpi_ensemble_.ensemble_communicator());
725 VectorTools::integrate_difference(
726 offline_data_.dof_handler(),
728 Functions::ZeroFunction<dim, Number>(),
731 VectorTools::L2_norm);
733 l2_norm_analytic = Number(std::sqrt(
734 Utilities::MPI::sum(
std::pow(difference_per_cell.l2_norm(), 2),
735 mpi_ensemble_.ensemble_communicator())));
740 U.extract_component(error_component, index);
742 offline_data_.affine_constraints().distribute(error_component);
743 error_component.update_ghost_values();
744 error_component -= analytic_component;
746 const Number linf_norm_error = Utilities::MPI::max(
747 error_component.linfty_norm(), mpi_ensemble_.ensemble_communicator());
749 VectorTools::integrate_difference(offline_data_.dof_handler(),
751 Functions::ZeroFunction<dim, Number>(),
754 VectorTools::L1_norm);
756 const Number l1_norm_error = Utilities::MPI::sum(
757 difference_per_cell.l1_norm(), mpi_ensemble_.ensemble_communicator());
759 VectorTools::integrate_difference(offline_data_.dof_handler(),
761 Functions::ZeroFunction<dim, Number>(),
764 VectorTools::L2_norm);
766 const Number l2_norm_error = Number(std::sqrt(
767 Utilities::MPI::sum(
std::pow(difference_per_cell.l2_norm(), 2),
768 mpi_ensemble_.ensemble_communicator())));
770 if (error_normalize_) {
771 linf_norm += linf_norm_error / linf_norm_analytic;
772 l1_norm += l1_norm_error / l1_norm_analytic;
773 l2_norm += l2_norm_error / l2_norm_analytic;
775 linf_norm += linf_norm_error;
776 l1_norm += l1_norm_error;
777 l2_norm += l2_norm_error;
781 if (mpi_ensemble_.ensemble_rank() != 0)
789 if (mpi_ensemble_.n_ensembles() > 1) {
790 linf_norm = Utilities::MPI::sum(
791 linf_norm, mpi_ensemble_.ensemble_leader_communicator());
792 l1_norm = Utilities::MPI::sum(
793 l1_norm, mpi_ensemble_.ensemble_leader_communicator());
794 l2_norm = Utilities::MPI::sum(
795 l2_norm, mpi_ensemble_.ensemble_leader_communicator());
798 if (mpi_ensemble_.world_rank() != 0)
801 logfile_ << std::endl <<
"Computed errors:" << std::endl << std::endl;
802 logfile_ << std::setprecision(16);
804 std::string description =
805 error_normalize_ ?
"Normalized consolidated" :
"Consolidated";
807 logfile_ << description +
" Linf, L1, and L2 errors at final time \n";
808 logfile_ << std::setprecision(16);
809 logfile_ <<
"#dofs = " << n_global_dofs_ << std::endl;
810 logfile_ <<
"t = " << t << std::endl;
811 logfile_ <<
"Linf = " << linf_norm << std::endl;
812 logfile_ <<
"L1 = " << l1_norm << std::endl;
813 logfile_ <<
"L2 = " << l2_norm << std::endl;
815 std::cout << description +
" Linf, L1, and L2 errors at final time \n";
816 std::cout << std::setprecision(16);
817 std::cout <<
"#dofs = " << n_global_dofs_ << std::endl;
818 std::cout <<
"t = " << t << std::endl;
819 std::cout <<
"Linf = " << linf_norm << std::endl;
820 std::cout <<
"L1 = " << l1_norm << std::endl;
821 std::cout <<
"L2 = " << l2_norm << std::endl;
825 template <
typename Description,
int dim,
typename Number>
827 const std::string &name,
829 const unsigned int cycle)
832 std::cout <<
"TimeLoop<dim, Number>::output(t = " << t <<
")" << std::endl;
835 const bool do_full_output =
836 (cycle % timer_output_full_multiplier_ == 0) && enable_output_full_;
837 const bool do_levelsets =
838 (cycle % timer_output_levelsets_multiplier_ == 0) &&
839 enable_output_levelsets_;
840 const bool do_checkpointing =
841 (cycle % timer_checkpoint_multiplier_ == 0) && enable_checkpointing_;
844 if (!(do_full_output || do_levelsets || do_checkpointing))
847 hyperbolic_module_.prepare_state_vector(state_vector, t);
848 if (!ParabolicSystem::is_identity)
849 parabolic_module_.prepare_state_vector(state_vector, t);
852 if (do_full_output || do_levelsets) {
853 Scope scope(computing_timer_,
"time step [X] - perform vtu output");
854 print_info(
"scheduling output");
856 postprocessor_.compute(state_vector);
863 postprocessor_.reset_bounds();
865 vtu_output_.schedule_output(
866 state_vector, name, t, cycle, do_full_output, do_levelsets);
870 if (do_checkpointing) {
871 Scope scope(computing_timer_,
"time step [X] - perform checkpointing");
872 print_info(
"scheduling checkpointing");
873 write_checkpoint(state_vector, base_name_ensemble_, t, cycle);
883 template <
typename Description,
int dim,
typename Number>
887 if (mpi_ensemble_.world_rank() != 0)
896 stream << std::endl <<
"Run time parameters:" << std::endl << std::endl;
897 ParameterAcceptor::prm.print_parameters(
898 stream, ParameterHandler::OutputStyle::ShortPRM);
903 std::ofstream output(base_name_ +
"-parameters.prm");
904 ParameterAcceptor::prm.print_parameters(output, ParameterHandler::ShortPRM);
908 template <
typename Description,
int dim,
typename Number>
919 std::vector<double> values = {
920 (double)offline_data_.n_export_indices(),
921 (double)offline_data_.n_locally_internal(),
922 (double)offline_data_.n_locally_owned(),
923 (double)offline_data_.n_locally_relevant(),
924 (double)offline_data_.n_export_indices() /
925 (double)offline_data_.n_locally_relevant(),
926 (double)offline_data_.n_locally_internal() /
927 (double)offline_data_.n_locally_relevant(),
928 (double)offline_data_.n_locally_owned() /
929 (double)offline_data_.n_locally_relevant()};
933 Utilities::MPI::min_max_avg(values, mpi_ensemble_.world_communicator());
935 if (mpi_ensemble_.world_rank() != 0)
938 std::ostringstream output;
941 dealii::Utilities::needed_digits(mpi_ensemble_.n_world_ranks());
943 const auto print_snippet = [&output, n](
const std::string &name,
944 const auto &values) {
945 output << name <<
": ";
947 output << std::setw(9) << (
unsigned int)values.min
948 <<
" [p" << std::setw(n) << values.min_index <<
"] "
949 << std::setw(9) << (
unsigned int)values.avg <<
" "
950 << std::setw(9) << (
unsigned int)values.max
951 <<
" [p" << std::setw(n) << values.max_index <<
"]";
955 const auto print_percentages = [&output, n](
const auto &percentages) {
956 output << std::endl <<
" ";
957 output <<
" (" << std::setw(3) << std::setprecision(2)
958 << percentages.min * 100 <<
"% )"
959 <<
" [p" << std::setw(n) << percentages.min_index <<
"] "
960 <<
" (" << std::setw(3) << std::setprecision(2)
961 << percentages.avg * 100 <<
"% )"
963 <<
" (" << std::setw(3) << std::setprecision(2)
964 << percentages.max * 100 <<
"% )"
965 <<
" [p" << std::setw(n) << percentages.max_index <<
"]";
968 output << std::endl << std::endl <<
"Partition: ";
969 print_snippet(
"exp", data[0]);
970 print_percentages(data[4]);
972 output << std::endl <<
" ";
973 print_snippet(
"int", data[1]);
974 print_percentages(data[5]);
976 output << std::endl <<
" ";
977 print_snippet(
"own", data[2]);
978 print_percentages(data[6]);
980 output << std::endl <<
" ";
981 print_snippet(
"rel", data[3]);
983 stream << output.str() << std::endl;
987 template <
typename Description,
int dim,
typename Number>
989 std::ostream &stream)
991 Utilities::System::MemoryStats stats;
992 Utilities::System::get_memory_stats(stats);
994 Utilities::MPI::MinMaxAvg data = Utilities::MPI::min_max_avg(
995 stats.VmRSS / 1024., mpi_ensemble_.world_communicator());
997 if (mpi_ensemble_.world_rank() != 0)
1000 std::ostringstream output;
1003 dealii::Utilities::needed_digits(mpi_ensemble_.n_world_ranks());
1005 output <<
"\nMemory: [MiB]"
1006 << std::setw(8) << data.min
1007 <<
" [p" << std::setw(n) << data.min_index <<
"] "
1008 << std::setw(8) << data.avg <<
" "
1009 << std::setw(8) << data.max
1010 <<
" [p" << std::setw(n) << data.max_index <<
"]";
1012 stream << output.str() << std::endl;
1016 template <
typename Description,
int dim,
typename Number>
1019 std::vector<std::ostringstream> output(computing_timer_.size());
1021 const auto equalize = [&]() {
1023 std::max_element(output.begin(),
1025 [](
const auto &left,
const auto &right) {
1026 return left.str().length() < right.str().length();
1028 const auto length = ptr->str().length();
1029 for (
auto &it : output)
1030 it << std::string(length - it.str().length() + 1,
' ');
1033 const auto print_wall_time = [&](
auto &timer,
auto &stream) {
1034 const auto wall_time = Utilities::MPI::min_max_avg(
1035 timer.wall_time(), mpi_ensemble_.world_communicator());
1037 constexpr auto eps = std::numeric_limits<double>::epsilon();
1042 const auto skew_negative = std::max(
1043 100. * (wall_time.min - wall_time.avg) / wall_time.avg - eps, -99.9);
1044 const auto skew_positive = std::min(
1045 100. * (wall_time.max - wall_time.avg) / wall_time.avg + eps, 99.9);
1047 stream << std::setprecision(2) << std::fixed << std::setw(8)
1048 << wall_time.avg <<
"s [sk: " << std::setprecision(1)
1049 << std::setw(5) << std::fixed << skew_negative <<
"%/"
1050 << std::setw(4) << std::fixed << skew_positive <<
"%]";
1052 dealii::Utilities::needed_digits(mpi_ensemble_.n_world_ranks());
1053 stream <<
" [p" << std::setw(n) << wall_time.min_index <<
"/"
1054 << wall_time.max_index <<
"]";
1057 const auto cpu_time_statistics =
1058 Utilities::MPI::min_max_avg(computing_timer_[
"time loop"].cpu_time(),
1059 mpi_ensemble_.world_communicator());
1060 const double total_cpu_time = cpu_time_statistics.sum;
1062 const auto print_cpu_time =
1063 [&](
auto &timer,
auto &stream,
bool percentage) {
1064 const auto cpu_time = Utilities::MPI::min_max_avg(
1065 timer.cpu_time(), mpi_ensemble_.world_communicator());
1067 stream << std::setprecision(2) << std::fixed << std::setw(9)
1068 << cpu_time.sum <<
"s ";
1071 stream <<
"(" << std::setprecision(1) << std::setw(4)
1072 << 100. * cpu_time.sum / total_cpu_time <<
"%)";
1075 auto jt = output.begin();
1076 for (
auto &it : computing_timer_)
1077 *jt++ <<
" " << it.first;
1080 jt = output.begin();
1081 for (
auto &it : computing_timer_)
1082 print_wall_time(it.second, *jt++);
1085 jt = output.begin();
1086 bool compute_percentages =
false;
1087 for (
auto &it : computing_timer_) {
1088 print_cpu_time(it.second, *jt++, compute_percentages);
1089 if (it.first.starts_with(
"time loop"))
1090 compute_percentages =
true;
1094 if (mpi_ensemble_.world_rank() != 0)
1097 stream << std::endl <<
"Timer statistics:\n";
1098 for (
auto &it : output)
1099 stream << it.str() << std::endl;
1103 template <
typename Description,
int dim,
typename Number>
1105 unsigned int cycle, Number t, std::ostream &stream,
bool final_time)
1111 static struct Data {
1112 unsigned int cycle = 0;
1114 double cpu_time_sum = 0.;
1115 double cpu_time_avg = 0.;
1116 double cpu_time_min = 0.;
1117 double cpu_time_max = 0.;
1118 double wall_time = 0.;
1119 } previous, current;
1121 static double time_per_second_exp = 0.;
1128 current.cycle = cycle;
1131 const auto wall_time_statistics =
1132 Utilities::MPI::min_max_avg(computing_timer_[
"time loop"].wall_time(),
1133 mpi_ensemble_.world_communicator());
1134 current.wall_time = wall_time_statistics.max;
1136 const auto cpu_time_statistics =
1137 Utilities::MPI::min_max_avg(computing_timer_[
"time loop"].cpu_time(),
1138 mpi_ensemble_.world_communicator());
1139 current.cpu_time_sum = cpu_time_statistics.sum;
1140 current.cpu_time_avg = cpu_time_statistics.avg;
1141 current.cpu_time_min = cpu_time_statistics.min;
1142 current.cpu_time_max = cpu_time_statistics.max;
1150 double delta_cycles = current.cycle - previous.cycle;
1151 const double cycles_per_second =
1152 delta_cycles / (current.wall_time - previous.wall_time);
1154 const auto efficiency = time_integrator_.efficiency();
1155 const auto n_dofs =
static_cast<double>(n_global_dofs_);
1157 const double wall_m_dofs_per_sec =
1158 delta_cycles * n_dofs / 1.e6 /
1159 (current.wall_time - previous.wall_time) * efficiency;
1161 double cpu_m_dofs_per_sec = delta_cycles * n_dofs / 1.e6 /
1162 (current.cpu_time_sum - previous.cpu_time_sum) *
1165 if (terminal_show_rank_throughput_)
1166 cpu_m_dofs_per_sec *= MultithreadInfo::n_threads();
1169 double cpu_time_skew = (current.cpu_time_max - current.cpu_time_min -
1170 previous.cpu_time_max + previous.cpu_time_min) /
1173 cpu_time_skew = std::max(0., cpu_time_skew);
1175 const double cpu_time_skew_percentage =
1176 cpu_time_skew * delta_cycles /
1177 (current.cpu_time_avg - previous.cpu_time_avg);
1179 const double delta_time =
1180 (current.t - previous.t) / (current.cycle - previous.cycle);
1181 const double time_per_second =
1182 (current.t - previous.t) / (current.wall_time - previous.wall_time);
1186 std::ostringstream output;
1189 output << std::endl;
1191 output <<
"Throughput:\n "
1192 << (terminal_show_rank_throughput_?
"RANK: " :
"CPU : ")
1193 << std::setprecision(4) << std::fixed << cpu_m_dofs_per_sec
1195 << std::scientific << 1. / cpu_m_dofs_per_sec * 1.e-6
1196 <<
" s/Qdof/substep)" << std::endl;
1198 output <<
" [cpu time skew: "
1199 << std::setprecision(2) << std::scientific << cpu_time_skew
1201 << std::setprecision(1) << std::setw(4) << std::setfill(
' ') << std::fixed
1202 << 100. * cpu_time_skew_percentage
1203 <<
"%)]" << std::endl;
1206 << std::setprecision(4) << std::fixed << wall_m_dofs_per_sec
1208 << std::scientific << 1. / wall_m_dofs_per_sec * 1.e-6
1209 <<
" s/Qdof/substep) ("
1210 << std::setprecision(2) << std::fixed << cycles_per_second
1211 <<
" cycles/s)" << std::endl;
1213 const auto &scheme = time_integrator_.time_stepping_scheme();
1215 << Patterns::Tools::Convert<TimeSteppingScheme>::to_string(scheme)
1217 << std::setprecision(2) << std::fixed << hyperbolic_module_.cfl()
1219 << std::setprecision(0) << std::fixed << hyperbolic_module_.n_restarts()
1221 << std::setprecision(0) << std::fixed << parabolic_module_.n_restarts()
1223 << std::setprecision(0) << std::fixed << hyperbolic_module_.n_warnings()
1225 << std::setprecision(0) << std::fixed << parabolic_module_.n_warnings()
1227 << std::setprecision(0) << std::fixed << hyperbolic_module_.n_corrections()
1229 << std::setprecision(0) << std::fixed << parabolic_module_.n_corrections()
1230 <<
" corr) ]" << std::endl;
1232 if constexpr (!ParabolicSystem::is_identity)
1233 parabolic_module_.print_solver_statistics(output);
1235 output <<
" [ dt = "
1236 << std::scientific << std::setprecision(2) << delta_time
1239 <<
" dt/s) ]" << std::endl;
1243 time_per_second_exp = 0.8 * time_per_second_exp + 0.2 * time_per_second;
1244 auto eta =
static_cast<unsigned int>(std::max(t_final_ - t, Number(0.)) /
1245 time_per_second_exp);
1247 output <<
"\n ETA : ";
1249 const unsigned int days = eta / (24 * 3600);
1251 output << days <<
" d ";
1255 const unsigned int hours = eta / 3600;
1257 output << hours <<
" h ";
1261 const unsigned int minutes = eta / 60;
1262 output << minutes <<
" min";
1264 if (mpi_ensemble_.world_rank() != 0)
1267 stream << output.str() << std::endl;
1271 template <
typename Description,
int dim,
typename Number>
1274 if (mpi_ensemble_.world_rank() != 0)
1277 std::cout <<
"[INFO] " << header << std::endl;
1281 template <
typename Description,
int dim,
typename Number>
1284 const std::string &secondary,
1285 std::ostream &stream)
1287 if (mpi_ensemble_.world_rank() != 0)
1290 const int header_size = header.size();
1291 const auto padded_header =
1292 std::string(std::max(0, 34 - header_size) / 2,
' ') + header +
1293 std::string(std::max(0, 35 - header_size) / 2,
' ');
1295 const int secondary_size = secondary.size();
1296 const auto padded_secondary =
1297 std::string(std::max(0, 34 - secondary_size) / 2,
' ') + secondary +
1298 std::string(std::max(0, 35 - secondary_size) / 2,
' ');
1302 stream <<
" ####################################################\n";
1303 stream <<
" #########" << padded_header <<
"#########\n";
1304 stream <<
" #########" << padded_secondary <<
"#########\n";
1305 stream <<
" ####################################################\n";
1306 stream << std::endl;
1311 template <
typename Description,
int dim,
typename Number>
1315 unsigned int timer_cycle,
1316 bool write_to_logfile,
1319 static const std::string vectorization_name = [] {
1320 constexpr auto width = VectorizedArray<Number>::size();
1326 result = std::to_string(width * 8 *
sizeof(Number)) +
" bit packed ";
1328 if constexpr (std::is_same_v<Number, double>)
1329 return result +
"double";
1330 else if constexpr (std::is_same_v<Number, float>)
1331 return result +
"float";
1336 std::ostringstream output;
1338 std::ostringstream primary;
1340 primary <<
"FINAL (cycle " << Utilities::int_to_string(cycle, 6) <<
")";
1342 primary <<
"Cycle " << Utilities::int_to_string(cycle, 6)
1343 <<
" (" << std::fixed << std::setprecision(1)
1344 << t / t_final_ * 100 <<
"%)";
1347 std::ostringstream secondary;
1348 secondary <<
"at time t = " << std::setprecision(8) << std::fixed << t;
1350 print_head(primary.str(), secondary.str(), output);
1352 output <<
"Information: (HYP) " << hyperbolic_system_.get().problem_name;
1353 if constexpr (!ParabolicSystem::is_identity) {
1354 output <<
"\n (PAR) " << parabolic_system_.get().problem_name;
1356 output <<
"\n [" << base_name_ <<
"] ";
1357 if (mpi_ensemble_.n_ensembles() > 1) {
1358 output << mpi_ensemble_.n_ensembles() <<
" ensembles ";
1361 << n_global_dofs_ <<
" Qdofs on "
1362 << mpi_ensemble_.n_world_ranks() <<
" ranks / "
1364 << MultithreadInfo::n_threads() <<
" threads <"
1366 <<
"[openmp disabled] <"
1368 << vectorization_name
1369 <<
">\n Last output cycle "
1371 <<
" at t = " << timer_granularity_ * (timer_cycle - 1)
1372 <<
" (terminal update interval " << terminal_update_interval_
1375 print_memory_statistics(output);
1376 print_timers(output);
1377 print_throughput(cycle, t, output, final_time);
1379 if (mpi_ensemble_.world_rank() == 0) {
1381 std::cout <<
"\033[2J\033[H";
1383 std::cout << output.str() << std::flush;
1385 if (write_to_logfile) {
1386 logfile_ <<
"\n" << output.str() << std::flush;
void write_checkpoint(const StateVector &state_vector, const std::string &base_name, const Number &t, const unsigned int &output_cycle)
Vectors::ScalarVector< Number > ScalarVector
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)
typename View::StateVector StateVector
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_cycle_statistics(unsigned int cycle, Number t, unsigned int output_cycle, bool write_to_logfile=false, 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)