diff --git a/README.md b/README.md index 34d4a064..36925602 100644 --- a/README.md +++ b/README.md @@ -184,13 +184,13 @@ The following options are available: (default value: 100) * newton\_tolerance: tolerance of the Newton solver (default value: 1e-6) * jfnk: use Jacobian-Free Newton Krylov method (default value: false) -* experiment: (optional) +* experiment (optional): * read\_in\_experimental\_data: whether to read in experimental data (default: false) * if reading in experimental data: * file: format of the file names. The format is pretty arbitrary, the keywords \#frame and \#camera are replaced by the frame and the camera number. The format of the file itself should be csv with a header line. (required) - * format: The format of the experimental data, either `point_cloud`, with (x,y,z,value) per line, or `ray`, with (pt0_x,pt0_y,pt0_z,pt1_x,pt1_y,pt1_z,value) per line, where the ray starts at pt0 and passes through pt1. (required) + * format: The format of the experimental data, either `point_cloud`, with (x,y,z,value) per line, or `ray`, with (pt0\_x,pt0\_y,pt0\_z,pt1\_x,pt1\_y,pt1\_z,value) per line, where the ray starts at pt0 and passes through pt1. (required) * first\_frame: number associated to the first frame (default value: 0) * last\_frame: number associated to the last frame (required) * first\_camera\_id: number associated to the first camera (required) @@ -203,19 +203,19 @@ The following options are available: given by a standard deviation (under the simplifying assumption that the error is normally distributed and independent for each data point) (default value: 0.0). * output\_experiment\_on\_mesh: Whether to output the experimental data projected onto the simulation mesh at each experiment time stamp (default: true). -* ensemble: (optional) +* ensemble (optional): * ensemble\_simulation: whether to perform an ensemble of simulations (default value: false) * ensemble\_size: the number of ensemble members for the ensemble Kalman filter (EnKF) (default value: 5) * initial\_temperature\_stddev: the standard deviation for the initial temperature of the material (default value: 0.0) * new\_material\_temperature\_stddev: the standard deviation for the temperature of material added during the process (default value: 0.0) * beam\_0\_max\_power\_stddev: the standard deviation for the max power for beam 0 (if it exists) (default value: 0.0) * beam\_0\_absorption\_efficiency\_stddev: the standard deviation for the absorption efficiency for beam 0 (if it exists) (default value: 0.0) -* data\_assimilation: (optional) +* data\_assimilation (optional): * assimilate\_data: whether to perform data assimilation (default value: false) * localization\_cutoff\_function: the function used to decrease the sample covariance as the relevant points become farther away: gaspari\_cohn, step\_function, none (default: none) * localization\_cutoff\_distance: the distance at which sample covariance entries are set to zero (default: infinity) * augment\_with\_beam\_0\_absorption: whether to augment the state vector with the beam 0 absorption efficiency (default: false) - * augment\_with\_beam\_0\_max_power: whether to augment the state vector with the beam 0 max power (default: false) + * augment\_with\_beam\_0\_max\_power: whether to augment the state vector with the beam 0 max power (default: false) * solver: * max\_number\_of\_temp\_vectors: maximum number of temporary vectors for the GMRES solve (optional) * max\_iterations: maximum number of iterations for the GMRES solve (optional) @@ -223,7 +223,13 @@ The following options are available: * profiling (optional): * timer: output timing information (default value: false) * caliper: configuration string for Caliper (optional) -* verbose_output: true or false (default value: false) +* checkpoint (optional): + * time\_steps\_between\_checkpoint: number of time steps after which + checkpointing is performed (required) + * filename\_prefix: prefix of the checkpoint files (required) +* restart (optional): + * filename\_prefix: prefix of the restart files (required) +* verbose\_output: true or false (default value: false) diff --git a/application/adamantine.hh b/application/adamantine.hh index 2592b8fa..98cafd51 100644 --- a/application/adamantine.hh +++ b/application/adamantine.hh @@ -37,9 +37,14 @@ #include #include +#include +#include #include +#include +#include #include +#include #include #include @@ -471,7 +476,7 @@ void refine_and_transfer( // Update the AffineConstraints and resize the solution thermal_physics->setup_dofs(); - thermal_physics->initialize_dof_vector(solution); + thermal_physics->initialize_dof_vector(0., solution); // Update MaterialProperty DoFHandler and resize the state vectors material_properties.reinit_dofs(); @@ -943,20 +948,63 @@ run(MPI_Comm const &communicator, boost::property_tree::ptree const &database, : mechanical_physics->get_dof_handler()); } + // Get the checkpoint and restart subtrees + boost::optional + checkpoint_optional_database = database.get_child_optional("checkpoint"); + boost::optional + restart_optional_database = database.get_child_optional("restart"); + + unsigned int time_steps_checkpoint = std::numeric_limits::max(); + std::string checkpoint_filename; + if (checkpoint_optional_database) + { + auto checkpoint_database = checkpoint_optional_database.get(); + // PropertyTreeInput checkpoint.time_steps_between_checkpoint + time_steps_checkpoint = + checkpoint_database.get("time_steps_between_checkpoint"); + // PropertyTreeInput checkpoint.filename_prefix + checkpoint_filename = + checkpoint_database.get("filename_prefix"); + } + + bool restart = false; + std::string restart_filename; + if (restart_optional_database) + { + restart = true; + auto restart_database = restart_optional_database.get(); + // PropertyTreeInput restart.filename_prefix + restart_filename = restart_database.get("filename_prefix"); + } + dealii::LA::distributed::Vector temperature; dealii::LA::distributed::Vector displacement; if (use_thermal_physics) { - thermal_physics->setup_dofs(); - thermal_physics->update_material_deposition_orientation(); - thermal_physics->compute_inverse_mass_matrix(); - thermal_physics->initialize_dof_vector(initial_temperature, temperature); - thermal_physics->get_state_from_material_properties(); + if (restart == false) + { + thermal_physics->setup(); + thermal_physics->initialize_dof_vector(initial_temperature, temperature); + } + else + { +#ifdef ADAMANTINE_WITH_CALIPER + CALI_MARK_BEGIN("restart from file"); +#endif + thermal_physics->load_checkpoint(restart_filename, temperature); +#ifdef ADAMANTINE_WITH_CALIPER + CALI_MARK_END("restart from file"); +#endif + } } if (use_mechanical_physics) { + // We currently do not support restarting the mechanical simulation. + adamantine::ASSERT_THROW( + restart == false, + "Mechanical simulation cannot be restarted from a file"); if (use_thermal_physics) { // Thermo-mechanical simulation @@ -979,6 +1027,19 @@ run(MPI_Comm const &communicator, boost::property_tree::ptree const &database, unsigned int n_time_step = 0; double time = 0.; double activation_time_end = -1.; + unsigned int const rank = + dealii::Utilities::MPI::this_mpi_process(communicator); + if (restart == true) + { + if (rank == 0) + { + std::cout << "Restarting from file" << std::endl; + } + std::ifstream file{restart_filename + "_time.txt"}; + boost::archive::text_iarchive ia{file}; + ia >> time; + ia >> n_time_step; + } // PropertyTreeInput geometry.deposition_time double const activation_time = geometry_database.get("deposition_time", 0.); @@ -1025,7 +1086,6 @@ run(MPI_Comm const &communicator, boost::property_tree::ptree const &database, #endif if ((time + time_step) > duration) time_step = duration - time; - unsigned int rank = dealii::Utilities::MPI::this_mpi_process(communicator); // Refine the mesh after time_steps_refinement time steps or when time // is greater or equal than the next predicted time for refinement. This @@ -1074,7 +1134,7 @@ run(MPI_Comm const &communicator, boost::property_tree::ptree const &database, // activation_end. timers[adamantine::add_material_search].start(); auto elements_to_activate = adamantine::get_elements_to_activate( - thermal_physics->get_dof_handler(), material_deposition_boxes); + thermal_physics->get_dof_handler(), material_deposition_boxes); timers[adamantine::add_material_search].stop(); // For now assume that all deposited material has never been melted @@ -1154,6 +1214,25 @@ run(MPI_Comm const &communicator, boost::property_tree::ptree const &database, time_step += time_step; } + if (n_time_step % time_steps_checkpoint == 0) + { +#ifdef ADAMANTINE_WITH_CALIPER + CALI_MARK_BEGIN("save checkpoint"); +#endif + if (rank == 0) + { + std::cout << "Checkpoint reached" << std::endl; + } + thermal_physics->save_checkpoint(checkpoint_filename, temperature); + std::ofstream file{checkpoint_filename + "_time.txt"}; + boost::archive::text_oarchive oa{file}; + oa << time; + oa << n_time_step; +#ifdef ADAMANTINE_WITH_CALIPER + CALI_MARK_END("save checkpoint"); +#endif + } + // Output progress on screen if (rank == 0) { @@ -1163,7 +1242,7 @@ run(MPI_Comm const &communicator, boost::property_tree::ptree const &database, if (int_part > progress) { std::cout << int_part * 10 << '%' << " completed" << std::endl; - ++progress; + progress = static_cast(int_part); } } @@ -1353,8 +1432,8 @@ run_ensemble(MPI_Comm const &global_communicator, ensemble_database.get("ensemble_size", 5); // Distribute the processors among the ensemble members MPI_Comm local_communicator; - unsigned int local_ensemble_size = -1; - unsigned int first_local_member = -1; + unsigned int local_ensemble_size = std::numeric_limits::max(); + unsigned int first_local_member = std::numeric_limits::max(); int my_color = -1; split_global_communicator(global_communicator, global_ensemble_size, local_communicator, local_ensemble_size, @@ -1462,6 +1541,37 @@ run_ensemble(MPI_Comm const &global_communicator, local_communicator, my_color, data_assimilation_database); + // Get the checkpoint and restart subtrees + boost::optional + checkpoint_optional_database = database.get_child_optional("checkpoint"); + boost::optional + restart_optional_database = database.get_child_optional("restart"); + + unsigned int const global_rank = + dealii::Utilities::MPI::this_mpi_process(global_communicator); + unsigned int time_steps_checkpoint = std::numeric_limits::max(); + std::string checkpoint_filename; + if (checkpoint_optional_database) + { + auto checkpoint_database = checkpoint_optional_database.get(); + // PropertyTreeInput checkpoint.time_steps_between_checkpoint + time_steps_checkpoint = + checkpoint_database.get("time_steps_between_checkpoint"); + // PropertyTreeInput checkpoint.filename_prefix + checkpoint_filename = + checkpoint_database.get("filename_prefix") + '_' + + std::to_string(global_rank); + } + + bool restart = false; + std::string restart_filename; + if (restart_optional_database) + { + restart = true; + auto restart_database = restart_optional_database.get(); + // PropertyTreeInput restart.filename_prefix + restart_filename = restart_database.get("filename_prefix"); + } for (unsigned int member = 0; member < local_ensemble_size; ++member) { // Resize the augmented ensemble block vector to have two blocks @@ -1525,18 +1635,27 @@ run_ensemble(MPI_Comm const &global_communicator, heat_sources_ensemble[member] = thermal_physics_ensemble[member]->get_heat_sources(); - thermal_physics_ensemble[member]->setup_dofs(); - thermal_physics_ensemble[member]->update_material_deposition_orientation(); - thermal_physics_ensemble[member]->compute_inverse_mass_matrix(); - - thermal_physics_ensemble[member]->initialize_dof_vector( - initial_temperature[member], - solution_augmented_ensemble[member].block(base_state)); - + if (restart == false) + { + thermal_physics_ensemble[member]->setup(); + thermal_physics_ensemble[member]->initialize_dof_vector( + initial_temperature[member], + solution_augmented_ensemble[member].block(base_state)); + } + else + { +#ifdef ADAMANTINE_WITH_CALIPER + CALI_MARK_BEGIN("restart from file"); +#endif + thermal_physics_ensemble[member]->load_checkpoint( + restart_filename + '_' + std::to_string(member), + solution_augmented_ensemble[member].block(base_state)); +#ifdef ADAMANTINE_WITH_CALIPER + CALI_MARK_END("restart from file"); +#endif + } solution_augmented_ensemble[member].collect_sizes(); - thermal_physics_ensemble[member]->get_state_from_material_properties(); - // For now we only output temperature post_processor_database.put("thermal_output", true); post_processor_ensemble.push_back( @@ -1558,11 +1677,10 @@ run_ensemble(MPI_Comm const &global_communicator, thermal_physics_ensemble[0]->get_dof_handler()); // ----- Read the experimental data ----- - unsigned int const global_rank = - dealii::Utilities::MPI::this_mpi_process(global_communicator); std::vector> frame_time_stamps; std::unique_ptr> experimental_data; - unsigned int experimental_frame_index = -1; + unsigned int experimental_frame_index = + std::numeric_limits::max(); if (experiment_optional_database) { @@ -1621,6 +1739,18 @@ run_ensemble(MPI_Comm const &global_communicator, unsigned int n_time_step = 0; double time = 0.; double activation_time_end = -1.; + if (restart == true) + { + if (global_rank == 0) + { + std::cout << "Restarting from file" << std::endl; + } + std::ifstream file{restart_filename + "_time.txt"}; + boost::archive::text_iarchive ia{file}; + ia >> time; + ia >> n_time_step; + } + // PropertyTreeInput geometry.deposition_time double const activation_time = geometry_database.get("deposition_time", 0.); @@ -1997,6 +2127,32 @@ run_ensemble(MPI_Comm const &global_communicator, } } + // ----- Checkpoint the ensemble members ----- + if (n_time_step % time_steps_checkpoint == 0) + { +#ifdef ADAMANTINE_WITH_CALIPER + CALI_MARK_BEGIN("save checkpoint"); +#endif + if (global_rank == 0) + { + std::cout << "Checkpoint reached" << std::endl; + } + for (unsigned int member = 0; member < local_ensemble_size; ++member) + { + thermal_physics_ensemble[member]->save_checkpoint( + checkpoint_filename + '_' + + std::to_string(first_local_member + member), + solution_augmented_ensemble[member].block(base_state)); + } + std::ofstream file{checkpoint_filename + "_time.txt"}; + boost::archive::text_oarchive oa{file}; + oa << time; + oa << n_time_step; +#ifdef ADAMANTINE_WITH_CALIPER + CALI_MARK_END("save checkpoint"); +#endif + } + // ----- Output progress on screen ----- if (global_rank == 0) { diff --git a/source/MaterialProperty.hh b/source/MaterialProperty.hh index 8cc129ba..2b3ac19a 100644 --- a/source/MaterialProperty.hh +++ b/source/MaterialProperty.hh @@ -1,4 +1,4 @@ -/* Copyright (c) 2016 - 2023, the adamantine authors. +/* Copyright (c) 2016 - 2024, the adamantine authors. * * This file is subject to the Modified BSD License and may not be distributed * without copyright and license information. Please refer to the file LICENSE @@ -210,6 +210,12 @@ public: std::vector> const &_cell_it_to_mf_pos, dealii::DoFHandler const &dof_handler); + /** + * Set the ratio of the material states at the cell level. + */ + void set_cell_state( + std::vector> const &cell_state); + /** * Return the underlying the DoFHandler. */ @@ -289,6 +295,8 @@ private: /** * Ratio of each in MaterarialState in each cell. */ + // FIXME Change the order of the indices. Currently, the first index is the + // state and the second is the cell. Kokkos::View _state; /** * Thermal properties of the material that are dependent of the state of the diff --git a/source/MaterialProperty.templates.hh b/source/MaterialProperty.templates.hh index 67010d8a..7a121de5 100644 --- a/source/MaterialProperty.templates.hh +++ b/source/MaterialProperty.templates.hh @@ -1,4 +1,4 @@ -/* Copyright (c) 2016 - 2023, the adamantine authors. +/* Copyright (c) 2016 - 2024, the adamantine authors. * * This file is subject to the Modified BSD License and may not be distributed * without copyright and license information. Please refer to the file LICENSE @@ -695,6 +695,22 @@ void MaterialProperty::set_state_device( }); } +template +void MaterialProperty::set_cell_state( + std::vector> const &cell_state) +{ + auto state_host = + Kokkos::create_mirror_view(Kokkos::WithoutInitializing, _state); + for (unsigned int i = 0; i < cell_state.size(); ++i) + { + for (unsigned int j = 0; j < g_n_material_states; ++j) + { + state_host(j, i) = cell_state[i][j]; + } + } + Kokkos::deep_copy(_state, state_host); +} + template void MaterialProperty::set_initial_state() { diff --git a/source/ThermalPhysics.hh b/source/ThermalPhysics.hh index 31a5aad8..caca5678 100644 --- a/source/ThermalPhysics.hh +++ b/source/ThermalPhysics.hh @@ -1,4 +1,4 @@ -/* Copyright (c) 2016 - 2023, the adamantine authors. +/* Copyright (c) 2016 - 2024, the adamantine authors. * * This file is subject to the Modified BSD License and may not be distributed * without copyright and license information. Please refer to the file LICENSE @@ -42,6 +42,8 @@ public: Geometry &geometry, MaterialProperty &material_properties); + void setup() override; + void setup_dofs() override; void compute_inverse_mass_matrix() override; @@ -72,13 +74,6 @@ public: double get_delta_t_guess() const override; - void initialize_dof_vector( - dealii::LA::distributed::Vector &vector) - const override; - - /** - * Initialize the given vector. The value is assumed to be a temperature. - */ void initialize_dof_vector(double const value, dealii::LA::distributed::Vector @@ -88,7 +83,13 @@ public: void set_state_to_material_properties() override; - void update_material_deposition_orientation() override; + void load_checkpoint(std::string const &filename, + dealii::LA::distributed::Vector + &temperature) override; + + void save_checkpoint(std::string const &filename, + dealii::LA::distributed::Vector + &temperature) override; void set_material_deposition_orientation( std::vector const &deposition_cos, @@ -125,6 +126,12 @@ private: using LA_Vector = typename dealii::LA::distributed::Vector; + /** + * Update the depostion cosine and sine from the Physics object to the + * operator object. + */ + void update_material_deposition_orientation(); + /** * Compute the right-hand side and apply the TermalOperator. */ diff --git a/source/ThermalPhysics.templates.hh b/source/ThermalPhysics.templates.hh index d4166fe5..1e0636bc 100644 --- a/source/ThermalPhysics.templates.hh +++ b/source/ThermalPhysics.templates.hh @@ -1,4 +1,4 @@ -/* Copyright (c) 2016 - 2023, the adamantine authors. +/* Copyright (c) 2016 - 2024, the adamantine authors. * * This file is subject to the Modified BSD License and may not be distributed * without copyright and license information. Please refer to the file LICENSE @@ -20,6 +20,7 @@ #include #include #include +#include #include #include #include @@ -532,6 +533,16 @@ ThermalPhysics::ThermalPhysics( _current_source_height = temp_height; } +template +void ThermalPhysics::setup() +{ + setup_dofs(); + update_material_deposition_orientation(); + compute_inverse_mass_matrix(); + get_state_from_material_properties(); +} + template void ThermalPhysics:: unsigned int const n_dofs_per_cell = _dof_handler.get_fe().n_dofs_per_cell(); unsigned int const direction_data_size = 2; unsigned int const phase_history_data_size = 1; - unsigned int constexpr n_material_states = - static_cast(adamantine::MaterialState::SIZE); unsigned int const data_size_per_cell = n_dofs_per_cell + direction_data_size + phase_history_data_size + - n_material_states; + g_n_material_states; dealii::Vector cell_solution(n_dofs_per_cell); std::vector dummy_cell_data(data_size_per_cell, std::numeric_limits::infinity()); @@ -920,15 +929,6 @@ double ThermalPhysics:: return time; } -template -void ThermalPhysics:: - initialize_dof_vector( - dealii::LA::distributed::Vector &vector) const -{ - _thermal_operator->initialize_dof_vector(vector); -} - template void ThermalPhysics:: @@ -936,10 +936,13 @@ void ThermalPhysics:: double const value, dealii::LA::distributed::Vector &vector) const { - // Resize the vector + // Resize the vector and initialize it to zero _thermal_operator->initialize_dof_vector(vector); - init_dof_vector(value, vector); + if (value != 0.) + { + init_dof_vector(value, vector); + } } template :: return solution; } + +template +void ThermalPhysics:: + load_checkpoint( + std::string const &filename, + dealii::LA::distributed::Vector &temperature) +{ + // Deserialize the mesh + auto &triangulation = _geometry.get_triangulation(); + triangulation.load(filename); + + // Deserialize the states, the direction, and the fe indices. + unsigned int const direction_data_size = 2; + unsigned int const data_size_per_cell = + g_n_material_states + direction_data_size + 1; + std::vector> data_to_deserialize( + triangulation.n_active_cells(), std::vector(data_size_per_cell)); + dealii::parallel::distributed::CellDataTransfer< + dim, dim, std::vector>> + cell_data_trans(triangulation); + cell_data_trans.deserialize(data_to_deserialize); + _deposition_cos.clear(); + _deposition_sin.clear(); + + unsigned int cell_id = 0; + std::vector> cell_state; + // Knowing the value of g_n_material_states simplifies the filling of + // cell_state but we need to make sure that the value of g_n_material_states + // is what we expect. + static_assert(g_n_material_states == 3); + for (auto const &cell : _dof_handler.active_cell_iterators()) + { + if (cell->is_locally_owned()) + { + // Get the state + cell_state.push_back({data_to_deserialize[cell_id][0], + data_to_deserialize[cell_id][1], + data_to_deserialize[cell_id][2]}); + + // Set the fe index + auto fe_index = static_cast( + data_to_deserialize[cell_id] + [g_n_material_states + direction_data_size]); + cell->set_active_fe_index(fe_index); + + // Get the direction + if (fe_index == 0) + { + _deposition_cos.push_back( + data_to_deserialize[cell_id][g_n_material_states]); + _deposition_sin.push_back( + data_to_deserialize[cell_id][g_n_material_states + 1]); + } + } + ++cell_id; + } + + setup_dofs(); + // Update MaterialProperty DoFHandler and resize the state vectors + _material_properties.reinit_dofs(); + // Update the state of each cell + _material_properties.set_cell_state(cell_state); + + // Finish the setup + _thermal_operator->set_material_deposition_orientation(_deposition_cos, + _deposition_sin); + compute_inverse_mass_matrix(); + get_state_from_material_properties(); + + // Deserialize the temperature + dealii::parallel::distributed::SolutionTransfer< + dim, dealii::LA::distributed::Vector> + solution_transfer(_dof_handler); + initialize_dof_vector(0., temperature); + if constexpr (std::is_same_v) + { + solution_transfer.deserialize(temperature); + } + else + { + dealii::LA::distributed::Vector + temperature_host(temperature.get_partitioner()); + solution_transfer.deserialize(temperature_host); + temperature.import_elements(temperature_host, + dealii::VectorOperation::insert); + } +} + +template +void ThermalPhysics:: + save_checkpoint( + std::string const &filename, + dealii::LA::distributed::Vector &temperature) +{ + // Prepare the states and the fe indices for serialization. + unsigned int const direction_data_size = 2; + unsigned int const data_size_per_cell = + g_n_material_states + direction_data_size + 1; + unsigned int locally_owned_cell_id = 0; + unsigned int activated_cell_id = 0; + unsigned int cell_id = 0; + auto &triangulation = _geometry.get_triangulation(); + std::vector> data_to_serialize( + triangulation.n_active_cells(), std::vector(data_size_per_cell)); + std::vector element_data(data_size_per_cell, 0.); + auto state_host = Kokkos::create_mirror_view_and_copy( + Kokkos::HostSpace{}, _material_properties.get_state()); + for (auto const &cell : _dof_handler.active_cell_iterators()) + { + if (cell->is_locally_owned()) + { + // Store the state + for (unsigned int i = 0; i < g_n_material_states; ++i) + { + data_to_serialize[cell_id][i] = state_host(i, locally_owned_cell_id); + } + + auto fe_index = cell->active_fe_index(); + // Store the direction + if (fe_index == 0) + { + data_to_serialize[cell_id][g_n_material_states] = + _deposition_cos[activated_cell_id]; + data_to_serialize[cell_id][g_n_material_states + 1] = + _deposition_sin[activated_cell_id]; + ++activated_cell_id; + } + else + { + // If there is no material, there is no deposition direction -> use an + // obviously wrong value. + data_to_serialize[cell_id][g_n_material_states] = 10.; + data_to_serialize[cell_id][g_n_material_states + 1] = 10.; + } + + // Store the FE index + data_to_serialize[cell_id][g_n_material_states + direction_data_size] = + fe_index; + + ++locally_owned_cell_id; + } + ++cell_id; + } + dealii::parallel::distributed::CellDataTransfer< + dim, dim, std::vector>> + cell_data_trans(triangulation); + cell_data_trans.prepare_for_serialization(data_to_serialize); + + // Prepare the temperature for serialization. We need to use a ghosted + // vector. + dealii::LA::distributed::Vector + ghosted_temperature( + temperature.locally_owned_elements(), + dealii::DoFTools::extract_locally_relevant_dofs(_dof_handler), + temperature.get_mpi_communicator()); + if constexpr (std::is_same_v) + { + ghosted_temperature = temperature; + } + else + { + dealii::LA::distributed::Vector + temperature_host(temperature.get_partitioner()); + temperature_host.import_elements(temperature, + dealii::VectorOperation::insert); + ghosted_temperature = temperature_host; + } + ghosted_temperature.update_ghost_values(); + dealii::parallel::distributed::SolutionTransfer< + dim, dealii::LA::distributed::Vector> + solution_transfer(_dof_handler); + solution_transfer.prepare_for_serialization(ghosted_temperature); + + // Serialize the mesh and the rest of the data. + triangulation.save(filename); +} } // namespace adamantine #endif diff --git a/source/ThermalPhysicsInterface.hh b/source/ThermalPhysicsInterface.hh index 90c30ce3..3c005346 100644 --- a/source/ThermalPhysicsInterface.hh +++ b/source/ThermalPhysicsInterface.hh @@ -1,4 +1,4 @@ -/* Copyright (c) 2016 - 2023, the adamantine authors. +/* Copyright (c) 2016 - 2024, the adamantine authors. * * This file is subject to the Modified BSD License and may not be distributed * without copyright and license information. Please refer to the file LICENSE @@ -36,6 +36,11 @@ public: virtual ~ThermalPhysicsInterface() = default; + /** + * Set up and initialize the data structure. + */ + virtual void setup() = 0; + /** * Associate the AffineConstraints and the MatrixFree objects to the * underlying Triangulation. @@ -84,13 +89,6 @@ public: */ virtual double get_delta_t_guess() const = 0; - /** - * Initialize the given vector. - */ - virtual void initialize_dof_vector( - dealii::LA::distributed::Vector &vector) - const = 0; - /** * Initialize the given vector with the given value. */ @@ -112,10 +110,20 @@ public: virtual void set_state_to_material_properties() = 0; /** - * Update the depostion cosine and sine from the Physics object to the - * operator object. + * Load the state of the simulation from files. */ - virtual void update_material_deposition_orientation() = 0; + virtual void + load_checkpoint(std::string const &filename, + dealii::LA::distributed::Vector + &temperature) = 0; + + /** + * Write the current state of the simulation on the filesystem. + */ + virtual void + save_checkpoint(std::string const &filename, + dealii::LA::distributed::Vector + &temperature) = 0; /** * Set the deposition cosine and sine and call diff --git a/tests/test_integration_3d.cc b/tests/test_integration_3d.cc index 6baa90f4..9ba3b037 100644 --- a/tests/test_integration_3d.cc +++ b/tests/test_integration_3d.cc @@ -94,3 +94,50 @@ BOOST_AUTO_TEST_CASE(integration_3D, *utf::tolerance(0.1)) } } } + +BOOST_AUTO_TEST_CASE(integration_3D_checkpoint_restart) +{ + MPI_Comm communicator = MPI_COMM_WORLD; + + std::vector timers; + initialize_timers(communicator, timers); + + // Read the input. + std::string const input_filename = "demo_316_short_anisotropic.info"; + adamantine::ASSERT_THROW(std::filesystem::exists(input_filename) == true, + "The file " + input_filename + " does not exist."); + boost::property_tree::ptree database; + boost::property_tree::info_parser::read_info(input_filename, database); + std::string const checkpoint_filename = "test_checkpoint"; + + // First run with checkpoint of the solution halfway through the solution + database.put("checkpoint.filename_prefix", checkpoint_filename); + database.put("checkpoint.time_steps_between_checkpoint", 80); + auto [temperature_1, displacement_1] = + run<3, dealii::MemorySpace::Host>(communicator, database, timers); + // Restart of the simulation + database.put("restart.filename_prefix", checkpoint_filename); + auto [temperature_2, displacement_2] = + run<3, dealii::MemorySpace::Host>(communicator, database, timers); + + // Compare the temperatures. When using more than one processor, the + // partitioning is different and so the distribution of the dofs are + // different. This makes it very difficult to compare the temperatures, so + // when using more than one processor we only compare the L2 norm. + if (dealii::Utilities::MPI::n_mpi_processes(communicator) == 1) + { + for (unsigned int i = 0; i < temperature_1.size(); ++i) + BOOST_TEST(temperature_1[i] == temperature_2[i]); + } + else + { + BOOST_TEST(temperature_1.l2_norm() == temperature_2.l2_norm()); + } + + // Remove the files created during the test + std::filesystem::remove(checkpoint_filename); + std::filesystem::remove(checkpoint_filename + "_fixed.data"); + std::filesystem::remove(checkpoint_filename + ".info"); + std::filesystem::remove(checkpoint_filename + "_variable.data"); + std::filesystem::remove(checkpoint_filename + "_time.txt"); +} diff --git a/tests/test_integration_3d_device.cc b/tests/test_integration_3d_device.cc index 360d1995..e2adfa6d 100644 --- a/tests/test_integration_3d_device.cc +++ b/tests/test_integration_3d_device.cc @@ -43,3 +43,50 @@ BOOST_AUTO_TEST_CASE(integration_3D_device, *utf::tolerance(0.1)) BOOST_TEST(temperature.local_element(i) == gold_value); } } + +BOOST_AUTO_TEST_CASE(integration_3D_checkpoint_restart_device) +{ + MPI_Comm communicator = MPI_COMM_WORLD; + + std::vector timers; + initialize_timers(communicator, timers); + + // Read the input. + std::string const input_filename = "demo_316_short_anisotropic.info"; + adamantine::ASSERT_THROW(std::filesystem::exists(input_filename) == true, + "The file " + input_filename + " does not exist."); + boost::property_tree::ptree database; + boost::property_tree::info_parser::read_info(input_filename, database); + std::string const checkpoint_filename = "test_checkpoint"; + + // First run with checkpoint of the solution halfway through the solution + database.put("checkpoint.filename_prefix", checkpoint_filename); + database.put("checkpoint.time_steps_between_checkpoint", 80); + auto [temperature_1, displacement_1] = + run<3, dealii::MemorySpace::Host>(communicator, database, timers); + // Restart of the simulation + database.put("restart.filename_prefix", checkpoint_filename); + auto [temperature_2, displacement_2] = + run<3, dealii::MemorySpace::Host>(communicator, database, timers); + + // Compare the temperatures. When using more than one processor, the + // partitioning is different and so the distribution of the dofs are + // different. This makes it very difficult to compare the temperatures, so + // when using more than one processor we only compare the L2 norm. + if (dealii::Utilities::MPI::n_mpi_processes(communicator) == 1) + { + for (unsigned int i = 0; i < temperature_1.size(); ++i) + BOOST_TEST(temperature_1[i] == temperature_2[i]); + } + else + { + BOOST_TEST(temperature_1.l2_norm() == temperature_2.l2_norm()); + } + + // Remove the files created during the test + std::filesystem::remove(checkpoint_filename); + std::filesystem::remove(checkpoint_filename + "_fixed.data"); + std::filesystem::remove(checkpoint_filename + ".info"); + std::filesystem::remove(checkpoint_filename + "_variable.data"); + std::filesystem::remove(checkpoint_filename + "_time.txt"); +} diff --git a/tests/test_material_deposition.cc b/tests/test_material_deposition.cc index 082f3906..4311b8b7 100644 --- a/tests/test_material_deposition.cc +++ b/tests/test_material_deposition.cc @@ -1,4 +1,4 @@ -/* Copyright (c) 2021 - 2022, the adamantine authors. +/* Copyright (c) 2021 - 2024, the adamantine authors. * * This file is subject to the Modified BSD License and may not be distributed * without copyright and license information. Please refer to the file LICENSE @@ -288,18 +288,15 @@ BOOST_AUTO_TEST_CASE(material_deposition) adamantine::ThermalPhysics> thermal_physics(communicator, database, geometry, material_properties); - thermal_physics.setup_dofs(); - thermal_physics.update_material_deposition_orientation(); - thermal_physics.compute_inverse_mass_matrix(); + thermal_physics.setup(); auto &dof_handler = thermal_physics.get_dof_handler(); auto [material_deposition_boxes, deposition_times, deposition_cos, deposition_sin] = adamantine::read_material_deposition(geometry_database); dealii::LA::distributed::Vector solution; + thermal_physics.initialize_dof_vector(0., solution); std::vector timers(adamantine::Timing::n_timers); - thermal_physics.initialize_dof_vector(solution); - thermal_physics.get_state_from_material_properties(); std::vector n_cells_ref = {610, 620, 630, 650, 650, 660, 670, 680, 720, 720}; double const time_step = 0.1; diff --git a/tests/test_mechanical_physics.cc b/tests/test_mechanical_physics.cc index 6a23fc54..9b6ea934 100644 --- a/tests/test_mechanical_physics.cc +++ b/tests/test_mechanical_physics.cc @@ -1,4 +1,4 @@ -/* Copyright (c) 2022 - 2023, the adamantine authors. +/* Copyright (c) 2022 - 2024, the adamantine authors. * * This file is subject to the Modified BSD License and may not be distributed * without copyright and license information. Please refer to the file LICENSE @@ -444,10 +444,7 @@ run_eshelby(std::vector> pts, unsigned int refinement_cycles) adamantine::ThermalPhysics> thermal_physics(communicator, database, geometry, material_properties); - thermal_physics.setup_dofs(); - thermal_physics.update_material_deposition_orientation(); - thermal_physics.compute_inverse_mass_matrix(); - thermal_physics.get_state_from_material_properties(); + thermal_physics.setup(); dealii::LinearAlgebra::distributed::Vector temperature; thermal_physics.initialize_dof_vector(100.0, temperature); diff --git a/tests/test_thermal_physics.hh b/tests/test_thermal_physics.hh index 98930f41..4eea10a7 100644 --- a/tests/test_thermal_physics.hh +++ b/tests/test_thermal_physics.hh @@ -1,4 +1,4 @@ -/* Copyright (c) 2016 - 2022, the adamantine authors. +/* Copyright (c) 2016 - 2024, the adamantine authors. * * This file is subject to the Modified BSD License and may not be distributed * without copyright and license information. Please refer to the file LICENSE @@ -127,13 +127,10 @@ void thermal_2d(boost::property_tree::ptree &database, double time_step) // Build ThermalPhysics adamantine::ThermalPhysics<2, 2, MemorySpaceType, dealii::QGauss<1>> physics( communicator, database, geometry, material_properties); - physics.setup_dofs(); - physics.update_material_deposition_orientation(); - physics.compute_inverse_mass_matrix(); - + physics.setup(); dealii::LA::distributed::Vector solution; - physics.initialize_dof_vector(solution); - physics.get_state_from_material_properties(); + physics.initialize_dof_vector(0., solution); + std::vector timers(adamantine::Timing::n_timers); double time = 0; while (time < 0.1) @@ -208,14 +205,11 @@ void thermal_2d_manufactured_solution() // Build ThermalPhysics adamantine::ThermalPhysics<2, 2, MemorySpaceType, dealii::QGauss<1>> physics( communicator, database, geometry, material_properties); - physics.setup_dofs(); - physics.update_material_deposition_orientation(); - physics.compute_inverse_mass_matrix(); - + physics.setup(); dealii::LA::distributed::Vector solution; + physics.initialize_dof_vector(0., solution); + std::vector timers(adamantine::Timing::n_timers); - physics.initialize_dof_vector(solution); - physics.get_state_from_material_properties(); double time = physics.evolve_one_time_step(0., 0.1, solution, timers); double const tolerance = 1e-5; @@ -260,13 +254,10 @@ void initial_temperature() // Build ThermalPhysics adamantine::ThermalPhysics<2, 2, MemorySpaceType, dealii::QGauss<1>> physics( communicator, database, geometry, material_properties); - physics.setup_dofs(); - physics.update_material_deposition_orientation(); - physics.compute_inverse_mass_matrix(); - + physics.setup(); dealii::LA::distributed::Vector solution; physics.initialize_dof_vector(1000., solution); - physics.get_state_from_material_properties(); + BOOST_TEST(solution.l1_norm() == 1000. * solution.size()); } @@ -325,15 +316,12 @@ void energy_conservation() // Build ThermalPhysics adamantine::ThermalPhysics<2, 2, MemorySpaceType, dealii::QGauss<1>> physics( communicator, database, geometry, material_properties); - physics.setup_dofs(); - physics.update_material_deposition_orientation(); - physics.compute_inverse_mass_matrix(); - + physics.setup(); dealii::LA::distributed::Vector solution; double constexpr initial_temperature = 10; double constexpr final_temperature = 10.5; physics.initialize_dof_vector(initial_temperature, solution); - physics.get_state_from_material_properties(); + std::vector timers(adamantine::Timing::n_timers); double time = 0; while (time < 100) @@ -439,14 +427,10 @@ void radiation_bcs() // Build ThermalPhysics adamantine::ThermalPhysics<2, 2, dealii::MemorySpace::Host, dealii::QGauss<1>> physics(communicator, database, geometry, material_properties); - physics.setup_dofs(); - physics.update_material_deposition_orientation(); - physics.compute_inverse_mass_matrix(); - + physics.setup(); dealii::LA::distributed::Vector solution; double constexpr initial_temperature = 10; physics.initialize_dof_vector(initial_temperature, solution); - physics.get_state_from_material_properties(); std::vector timers(adamantine::Timing::n_timers); double time = 0; while (time < 100) @@ -545,14 +529,10 @@ void convection_bcs() // Build ThermalPhysics adamantine::ThermalPhysics<2, 2, dealii::MemorySpace::Host, dealii::QGauss<1>> physics(communicator, database, geometry, material_properties); - physics.setup_dofs(); - physics.update_material_deposition_orientation(); - physics.compute_inverse_mass_matrix(); - + physics.setup(); dealii::LA::distributed::Vector solution; double constexpr initial_temperature = 10; physics.initialize_dof_vector(initial_temperature, solution); - physics.get_state_from_material_properties(); std::vector timers(adamantine::Timing::n_timers); double time = 0; while (time < 100) @@ -612,13 +592,9 @@ void reference_temperature() // Build ThermalPhysics adamantine::ThermalPhysics<2, 2, MemorySpaceType, dealii::QGauss<1>> physics( communicator, database, geometry, material_properties); - physics.setup_dofs(); - physics.update_material_deposition_orientation(); - physics.compute_inverse_mass_matrix(); - + physics.setup(); dealii::LA::distributed::Vector solution; physics.initialize_dof_vector(1000., solution); - physics.get_state_from_material_properties(); // Now check that the melting indicator works as expected std::vector reference_temperatures({1500.0, 300.0});