diff --git a/include/base/restart.h b/include/base/restart.h new file mode 100644 index 0000000..f22a131 --- /dev/null +++ b/include/base/restart.h @@ -0,0 +1,153 @@ +#pragma once + +#include + +#include +#include + +#include + +// for implementing a check-pointing and restart mechanism. +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace Utilities +{ + /** + * @brief Create files required for a restart consisting of the triangulation, vectors, and some metadata + * + * @tparam PartitionerPtr shred_ptr holding a fully distributed Partitioner with ghost elements + * + * @param triangulation The triangulation used to create the checkpoint + * @param dof_handler DoFHandeler associated to the global vectors + * @param partitioner A valid shared_ptr with partitioner used to create tmp vectors if necessary + * @param vectors A list of vectors, which are stored in the checkpoint + * @param name File name(s) for the restart files + * @param t Absolute time associated to the data vectors + * @param t Timestep associated to the data vectors + */ + template + void + create_restart_snapshot( + const dealii::parallel::distributed::Triangulation &triangulation, + const dealii::DoFHandler &dof_handler, + PartitionerPtr partitioner, + const std::vector &vectors, + const std::string &name, + const double t, + unsigned int timestep) + { + // To load a restart, we need read access to ghost entries. + // Thus, we might need temporary vectors to comply with the required + // partitioner layout (fully distributed vector with read access to ghost + // entries). We use lazy allocation, if required. + std::vector tmp_vectors(vectors.size()); + std::vector tmp_vectors_ptr; + + for (std::size_t i = 0; i < vectors.size(); ++i) + { + if (!vectors[i]->partitioners_are_globally_compatible( + *partitioner.get())) + { + // in case our vector carries anyway only local data, i.e., no ghost + // values at all, we need to create a temporary vector + tmp_vectors[i].reinit(partitioner); + tmp_vectors[i] = *(vectors[i]); + tmp_vectors_ptr.emplace_back(&tmp_vectors[i]); + } + else + { + // if the partitioner is compatible, we can just use the existing + // vector + tmp_vectors_ptr.emplace_back(vectors[i]); + } + } + + dealii::parallel::distributed::SolutionTransfer + solution_transfer(dof_handler); + + solution_transfer.prepare_for_serialization(tmp_vectors_ptr); + triangulation.save(name + ".mesh"); + + if (dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + { + std::ofstream file(name + ".metadata", std::ios::binary); + boost::archive::binary_oarchive oa(file); + oa << t; + oa << timestep; + } + } + + + /** + * @brief Load a restart snapshot previously stored using create_restart_snapshot + * + * @param dof_handler DoFHandler associated to the vectors + * @param vectors Vectors container to load data into + * @param name Name of the file bundle + * @param t Absolute time stored in the metadata + * @param timestep Timestep stored in the metadata + */ + template + void + load_restart_snapshot(const dealii::DoFHandler &dof_handler, + std::vector &vectors, + const std::string &name, + double &t, + unsigned int ×tep) + { + auto partitioner = std::make_shared( + dof_handler.locally_owned_dofs(), + DoFTools::extract_locally_relevant_dofs(dof_handler), + MPI_COMM_WORLD); + + // To load a restart, we need write access to ghost entries + // Thus, we might need temporary vectors to comply with the required + // partitioner layout (fully distributed vector with write access to ghost + // entries). We use lazy allocation, if required. + std::vector tmp_vectors(vectors.size()); + std::vector tmp_vectors_ptr; + + for (std::size_t i = 0; i < vectors.size(); ++i) + { + // With LA:d:V, we can grant write access using zero_out_ghost_values + // (this wouldn't work with distributed PETSc vectors or similar) + if (vectors[i]->has_ghost_elements()) + { + vectors[i]->zero_out_ghost_values(); + tmp_vectors_ptr.emplace_back(vectors[i]); + } + else + { + // in case our vector carries anyway only local data, i.e., no ghost + // values at all, we need to create a temporary vector + tmp_vectors[i].reinit(partitioner); + tmp_vectors_ptr.emplace_back(&tmp_vectors[i]); + } + } + + // triangulation.load(name + "-checkpoint.mesh"); + dealii::parallel::distributed::SolutionTransfer + solution_transfer(dof_handler); + solution_transfer.deserialize(tmp_vectors_ptr); + + for (unsigned int i = 0; i < tmp_vectors_ptr.size(); ++i) + { + // Copy over data from the temporary vectors, if necessary + if (!vectors[i]->has_ghost_elements()) + vectors[i]->copy_locally_owned_data_from(*(tmp_vectors_ptr[i])); + // ... and update all ghost values + vectors[i]->update_ghost_values(); + } + + // Last but not least, retrieve the time-stamp data + std::ifstream file(name + ".metadata", std::ios::binary); + boost::archive::binary_iarchive ia(file); + ia >> t; + ia >> timestep; + } +} // namespace Utilities + +DEAL_II_NAMESPACE_CLOSE diff --git a/include/base/time_handler.h b/include/base/time_handler.h index c7e29ef..d4a243b 100644 --- a/include/base/time_handler.h +++ b/include/base/time_handler.h @@ -45,6 +45,12 @@ class Time time_current += delta_t; ++timestep; } + void + set_time(double absolute_time, unsigned int nth_timestep) + { + time_current = absolute_time; + timestep = nth_timestep; + } private: unsigned int timestep; diff --git a/include/base/utilities.h b/include/base/utilities.h index b6c4ef3..8db5225 100644 --- a/include/base/utilities.h +++ b/include/base/utilities.h @@ -104,7 +104,42 @@ namespace Utilities << std::endl; } - + /** + * @brief Helper function to convert different representations of the same number + * (in double rpecision) to a common format by truncating trailing zeros. + * + * Examples: + * number = 0.0100 -> 0.01 + * number = 42.0100 -> 42.01 + * number = 42.00 -> 42 + * + * + * @param number the number to format + * @return std::string the string representation + */ + std::string + format_time_stamp_to_string(double number) + { + // We need to apply some cosmetics to t to make it usable + std::string str = std::to_string(number); + for (int i = str.size() - 1; i >= 1; i--) + { + if (str.at(i) == '0') + { + str.pop_back(); // Remove if last digit is '0'. + } + else if (str.at(i) == '.') + { + str.pop_back(); // Remove dot. + break; // Break after '.' is removed. + } + else + { + break; // Or break before a digit is removed. + } + } + return str; + } /** * @brief Rounds a given number up to a defined precision. diff --git a/include/parameter/parameter_handling.h b/include/parameter/parameter_handling.h index 8d0f638..8b9083b 100644 --- a/include/parameter/parameter_handling.h +++ b/include/parameter/parameter_handling.h @@ -311,8 +311,9 @@ namespace Parameters // Set the timestep size $ \varDelta t $ and the simulation end-time. struct Time { - double delta_t = 0.1; - double end_time = 1.; + double delta_t = 0.1; + double end_time = 1.; + double start_time = 0; void add_parameters(ParameterHandler &prm); @@ -323,6 +324,10 @@ namespace Parameters { prm.enter_subsection("Time"); { + prm.add_parameter("Start time", + start_time, + "Restart simulation or start from the beginning", + Patterns::Double(0.)); prm.add_parameter("End time", end_time, "End time", Patterns::Double()); prm.add_parameter("Time step size", diff --git a/include/solid_mechanics/mf_elasticity.h b/include/solid_mechanics/mf_elasticity.h index 6e0b010..2a93b47 100644 --- a/include/solid_mechanics/mf_elasticity.h +++ b/include/solid_mechanics/mf_elasticity.h @@ -62,6 +62,7 @@ static const unsigned int debug_level = 0; #include #include #include +#include #include #include #include @@ -79,8 +80,6 @@ using namespace dealii; // into it: namespace FSI { - using namespace dealii; - // The Solid class is the central class. template class Solid @@ -146,6 +145,13 @@ namespace FSI void output_results(const unsigned int result_number) const; + void + write_restart(); + + + void + load_restart(); + // Set up an Additional data object template void @@ -558,7 +564,8 @@ namespace FSI testcase = testcase_; make_grid(); system_setup(); - output_results(0); + if (parameters.start_time != 0) + output_results(0); time.increment(); // We then declare the incremental solution update $\varDelta @@ -573,6 +580,9 @@ namespace FSI mf_data_reference); precice_adapter->initialize(total_displacement); + if (parameters.start_time != 0) + load_restart(); + // At the beginning, we reset the solution update for this time step... while (precice_adapter->is_coupling_ongoing()) { @@ -614,6 +624,11 @@ namespace FSI time.increment(); } } + // write the restart after we finished + // we decrement the time in the write_restart function, as we call the + // function after we incremented the time for the next iteration + write_restart(); + // for post-processing, print average CG iterations over the whole run: timer_out << std::endl @@ -635,7 +650,18 @@ namespace FSI { Assert(testcase.get() != nullptr, ExcInternalError()); testcase->make_coarse_grid_and_bcs(triangulation); - triangulation.refine_global(parameters.n_global_refinement); + + if (parameters.start_time != 0) + { + std::string restart_file( + parameters.output_folder + "restart-t-" + + Utilities::format_time_stamp_to_string(parameters.start_time) + + ".mesh"); + pcout << "-- Looking for restart file \"" << restart_file << "\"\n"; + triangulation.load(restart_file); + } + else + triangulation.refine_global(parameters.n_global_refinement); if (testcase->body_force.get() == nullptr) testcase->body_force = @@ -2041,6 +2067,79 @@ namespace FSI return std::make_tuple(lin_it, lin_res, cond_number); } + + + template + void + Solid::write_restart() + { + // first, decrement the time, as we call the function after an (unnecessary) + // increment + double t_decrement = time.current() - time.get_delta_t(); + unsigned int timestep_decrement = time.get_timestep() - 1; + + auto t_string = Utilities::format_time_stamp_to_string(t_decrement); + + std::string restart_file(parameters.output_folder + "restart-t-" + + t_string); + pcout << "-- Creating restart files \"" + restart_file + + "\" for t = " + t_string + << " ( timestep " << timestep_decrement << " ) " + << "\n"; + std::vector in_vectors({&total_displacement, + &velocity, + &acceleration, + &old_displacement, + &velocity_old, + &acceleration_old}); + + // Make sure to pass updated vectors into the function + for (auto &in : in_vectors) + in->update_ghost_values(); + + Utilities::create_restart_snapshot( + triangulation, + dof_handler, + total_displacement.get_partitioner(), + in_vectors, + restart_file, + t_decrement, + timestep_decrement); + } + + + + template + void + Solid::load_restart() + { + std::vector in_vectors({&total_displacement, + &velocity, + &acceleration, + &old_displacement, + &velocity_old, + &acceleration_old}); + + std::string restart_file( + parameters.output_folder + "restart-t-" + + Utilities::format_time_stamp_to_string(parameters.start_time)); + + double loaded_time = 0; + unsigned int loaded_timestep = 0; + Utilities::load_restart_snapshot( + dof_handler, in_vectors, restart_file, loaded_time, loaded_timestep); + + pcout << "-- Restarting computation from files \"" + restart_file + + "\" at t = " + << loaded_time << " ( timestep " << loaded_timestep << " ) " + << "\n"; + time.set_time(loaded_time, loaded_timestep); + // we loaded the time we already computed, so let's move on to the next step + time.increment(); + } + + + template void Solid::output_results(const unsigned int result_number) const