From b25eb44fade01c847c32048f13fc8047b4542790 Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Sat, 31 Jul 2021 09:23:37 -0500 Subject: [PATCH] Step-69: Use SolutionTransfer instead of low-level boost archive --- examples/step-69/step-69.cc | 208 ++++++++++++++++++++---------------- 1 file changed, 113 insertions(+), 95 deletions(-) diff --git a/examples/step-69/step-69.cc b/examples/step-69/step-69.cc index 49edcf57a3..6ea0cadbf1 100644 --- a/examples/step-69/step-69.cc +++ b/examples/step-69/step-69.cc @@ -46,6 +46,7 @@ #include #include +#include #include #include @@ -162,6 +163,8 @@ namespace Step69 const QGauss quadrature; const QGauss face_quadrature; + unsigned int refinement; + private: TimerOutput &computing_timer; @@ -169,8 +172,6 @@ namespace Step69 double height; double disk_position; double disk_diameter; - - unsigned int refinement; }; // @sect4{The OfflineData class} @@ -258,7 +259,6 @@ namespace Step69 // $[\rho,\textbf{m},E]$. // // The purpose of the class members component_names, - // component_physical_units, // pressure, and speed_of_sound is evident from // their names. We also provide a function // compute_lambda_max(), that computes the wave speed @@ -293,8 +293,6 @@ namespace Step69 using flux_type = Tensor<1, problem_dimension, Tensor<1, dim>>; const static std::array component_names; - const static std::array - component_physical_units; static constexpr double gamma = 7. / 5.; @@ -486,11 +484,15 @@ namespace Step69 private: vector_type interpolate_initial_values(const double t = 0); + void checkpoint(const vector_type &U, + const std::string &name, + double t, + unsigned int cycle); + void output(const vector_type &U, const std::string &name, double t, - unsigned int cycle, - bool checkpoint = false); + unsigned int cycle); const MPI_Comm mpi_communicator; std::ostringstream timer_output; @@ -1600,43 +1602,23 @@ namespace Step69 } // We conclude this section by defining static arrays - // component_names and component_physical_units - // that contain strings describing the - // components of our state vector, as well as the physical units of - // these quantities that can be used to annotate the VTU output - // files we generate for visualization and postprocessing. We have - // template specializations for dimensions one, two and three, that - // are used later in DataOut for naming the corresponding - // components: + // component_names that contain strings describing the + // component names of our state vector. We have template specializations + // for dimensions one, two and three, that are used later in DataOut for + // naming the corresponding components: template <> const std::array ProblemDescription<1>::component_names{ {"rho", "m", "E"}}; - template <> - const std::array - ProblemDescription<1>::component_physical_units{ - {"kg/m", "kg*m/s/m", "J/m"}}; - - template <> const std::array ProblemDescription<2>::component_names{ {"rho", "m_1", "m_2", "E"}}; - template <> - const std::array - ProblemDescription<2>::component_physical_units{ - {"kg/m/m", "kg*m/s/m/m", "kg*m/s/m/m", "J/m/m"}}; - template <> const std::array ProblemDescription<3>::component_names{ {"rho", "m_1", "m_2", "m_3", "E"}}; - template <> - const std::array - ProblemDescription<3>::component_physical_units{ - {"kg/m/m/m", "kg*m/s/m/m/m", "kg*m/s/m/m/m", "kg*m/s/m/m/m", "J/m/m/m"}}; - // @sect4{Initial values} // The last preparatory step, before we discuss the implementation of the @@ -2487,11 +2469,28 @@ namespace Step69 pcout << "done" << std::endl; // Next we create the triangulation, assemble all matrices, set up - // scratch space, and initialize the DataOut object: + // scratch space, and initialize the DataOut object. All of these + // operations are pretty standard and discussed in detail in the + // Discretization and OfflineData classes. + // + // We have to make take care of a special case when resuming an + // interrupted computation though: In order to be able to read in the + // saved mesh and associated state vector we have to make sure to + // not refine the coarse mesh: { print_head(pcout, "create triangulation"); - discretization.setup(); + + if (resume) + { + discretization.refinement = 0; + discretization.setup(); + discretization.triangulation.load(base_name + "-checkpoint.mesh"); + } + else + { + discretization.setup(); + } pcout << "Number of active cells: " << discretization.triangulation.n_global_active_cells() @@ -2515,8 +2514,9 @@ namespace Step69 double t = 0.; unsigned int output_cycle = 0; - print_head(pcout, "interpolate initial values"); - vector_type U = interpolate_initial_values(); + vector_type U; + for (auto &it : U) + it.reinit(offline_data.partitioner); // @sect5{Resume} // @@ -2525,36 +2525,48 @@ namespace Step69 // resume==true we indicate that we have indeed an // interrupted computation and the program shall restart by reading in // an old state consisting of t, - // output_cycle, and U from a checkpoint - // file. These checkpoint files will be created in the - // output() routine discussed below. + // output_cycle, and the state vector U from + // checkpoint files. + // + // A this point we have already read in the stored refinement history + // of our parallel distributed mesh. What is missing are the actual + // state vector U, the time and output cycle. We use the + // SolutionTransfer class in combination with the + // distributed::Triangulation::load() / + // distributed::Triangulation::save() mechanism to read in the state + // vector. A separate boost::archive is used to retrieve + // t and output_cycle. The checkpoint files + // will be created in the output() routine discussed + // below. if (resume) { - print_head(pcout, "restore interrupted computation"); + print_head(pcout, "resume interrupted computation"); + + parallel::distributed:: + SolutionTransfer> + solution_transfer(offline_data.dof_handler); - const unsigned int i = - discretization.triangulation.locally_owned_subdomain(); + std::vector *> vectors; + std::transform(U.begin(), + U.end(), + std::back_inserter(vectors), + [](auto &it) { return ⁢ }); + solution_transfer.deserialize(vectors); - const std::string name = base_name + "-checkpoint-" + - Utilities::int_to_string(i, 4) + ".archive"; - std::ifstream file(name, std::ios::binary); + for (auto &it : U) + it.update_ghost_values(); - // We use a boost::archive to store and read in the - // contents the checkpointed state. + std::ifstream file(base_name + "-checkpoint.metadata", + std::ios::binary); boost::archive::binary_iarchive ia(file); ia >> t >> output_cycle; - - for (auto &it1 : U) - { - // it1 iterates over all components of the state - // vector U. We read in every entry of the - // component in sequence and update the ghost layer afterwards: - for (auto &it2 : it1) - ia >> it2; - it1.update_ghost_values(); - } + } + else + { + print_head(pcout, "interpolate initial values"); + U = interpolate_initial_values(); } // With either the initial state set up, or an interrupted state @@ -2594,7 +2606,8 @@ namespace Step69 if (t > output_cycle * output_granularity) { - output(U, base_name, t, output_cycle, true); + checkpoint(U, base_name, t, output_cycle); + output(U, base_name, t, output_cycle); ++output_cycle; } } @@ -2652,7 +2665,41 @@ namespace Step69 } // @sect5{Output and checkpointing} - // + + // We checkpoint the current state by doing the precise inverse + // operation to what we discussed for the resume + // logic: + + template + void MainLoop::checkpoint(const typename MainLoop::vector_type &U, + const std::string &name, + const double t, + const unsigned int cycle) + { + print_head(pcout, "checkpoint computation"); + + parallel::distributed:: + SolutionTransfer> + solution_transfer(offline_data.dof_handler); + + std::vector *> vectors; + std::transform(U.begin(), + U.end(), + std::back_inserter(vectors), + [](auto &it) { return ⁢ }); + + solution_transfer.prepare_for_serialization(vectors); + + discretization.triangulation.save(name + "-checkpoint.mesh"); + + if (dealii::Utilities::MPI::this_mpi_process(mpi_communicator) == 0) + { + std::ofstream file(name + "-checkpoint.metadata", std::ios::binary); + boost::archive::binary_oarchive oa(file); + oa << t << cycle; + } + } + // Writing out the final vtk files is quite an IO intensive task that can // stall the main loop for a while. In order to avoid this we use an asynchronous @@ -2670,11 +2717,9 @@ namespace Step69 void MainLoop::output(const typename MainLoop::vector_type &U, const std::string & name, const double t, - const unsigned int cycle, - const bool checkpoint) + const unsigned int cycle) { - pcout << "MainLoop::output(t = " << t - << ", checkpoint = " << checkpoint << ")" << std::endl; + pcout << "MainLoop::output(t = " << t << ")" << std::endl; // If the asynchronous writeback option is set we launch a background // thread performing all the slow IO to disc. In that case we have to @@ -2735,39 +2780,12 @@ namespace Step69 // the this pointer as well as most of the arguments of // the output function by value so that we have access to them inside // the lambda function. - const auto output_worker = [this, name, t, cycle, checkpoint, data_out]() { - if (checkpoint) - { - // We checkpoint the current state by doing the precise inverse - // operation to what we discussed for the resume - // logic: - - const unsigned int i = - discretization.triangulation.locally_owned_subdomain(); - std::string filename = - name + "-checkpoint-" + Utilities::int_to_string(i, 4) + ".archive"; - - std::ofstream file(filename, std::ios::binary | std::ios::trunc); - - boost::archive::binary_oarchive oa(file); - oa << t << cycle; - for (const auto &it1 : output_vector) - for (const auto &it2 : it1) - oa << it2; - } - - DataOutBase::VtkFlags output_flags; - output_flags.time = t; - output_flags.cycle = cycle; - output_flags.compression_level = DataOutBase::VtkFlags::best_speed; - for (unsigned int i = 0; - i < ProblemDescription::component_names.size(); - ++i) - output_flags - .physical_units[ProblemDescription::component_names[i]] = - ProblemDescription::component_physical_units[i]; - - data_out->set_flags(output_flags); + const auto output_worker = [this, name, t, cycle, data_out]() { + DataOutBase::VtkFlags flags(t, + cycle, + true, + DataOutBase::VtkFlags::best_speed); + data_out->set_flags(flags); data_out->write_vtu_with_pvtu_record( "", name + "-solution", cycle, mpi_communicator, 6); -- 2.39.5