From fbe5727f861cdf70f5799519cc7c6064f1c8e0fe Mon Sep 17 00:00:00 2001 From: mivaia Date: Fri, 2 Aug 2024 12:09:11 -0400 Subject: [PATCH] added set_solution_from_checkpoint function to clean the use of the restart with the initial condition. Scaling bug in average velocity --- include/solvers/navier_stokes_base.h | 8 +- source/solvers/gls_navier_stokes.cc | 31 ++-- source/solvers/navier_stokes_base.cc | 232 +++++++++++++-------------- 3 files changed, 140 insertions(+), 131 deletions(-) diff --git a/include/solvers/navier_stokes_base.h b/include/solvers/navier_stokes_base.h index a33cb34e8d..746e7de034 100644 --- a/include/solvers/navier_stokes_base.h +++ b/include/solvers/navier_stokes_base.h @@ -382,11 +382,17 @@ class NavierStokesBase : public PhysicsSolver refine_mesh_uniform(); /** - * @brief read_checkpoint + * @brief Restart a previous simulation from a checkpoint file. */ virtual void read_checkpoint(); + /** + * @brief Read the solution from a checkpoint file and set it as the current solution. + */ + virtual void + set_solution_from_checkpoint(std::string checkpoint_file_prefix); + /** * @brief set_nodal_values * Set the nodal values of velocity and pressure diff --git a/source/solvers/gls_navier_stokes.cc b/source/solvers/gls_navier_stokes.cc index b6952b4c63..9c9b95f407 100644 --- a/source/solvers/gls_navier_stokes.cc +++ b/source/solvers/gls_navier_stokes.cc @@ -1281,29 +1281,34 @@ GLSNavierStokesSolver::set_initial_condition_fd( else if (initial_condition_type == Parameters::InitialConditionType::average_velocity_profile) { - this->pcout << "*********************************************************" - << std::endl; - this->pcout - << " Initial condition using average velocity from checkpoint " - << std::endl; - this->pcout << "*********************************************************" - << std::endl; + announce_string(this->pcout, "Initial condition using average velocity profile"); + // this->pcout << "*********************************************************" + // << std::endl; + // this->pcout + // << " Initial condition using average velocity from checkpoint " + // << std::endl; + // this->pcout << "*********************************************************" + // << std::endl; + + + std::string checkpoint_file_name = + this->simulation_parameters.initial_condition->checkpoint_folder + this->simulation_parameters.initial_condition->checkpoint_file_name; // The average velocity profile needs to come from the results of a // simulation. The checkpoint/restart mechanism is used to obtain all the // simulations information necessary to restart the simulation. The same // mesh needs to be used, but new physics can be added to the simulation. - this->read_checkpoint(); + this->set_solution_from_checkpoint(checkpoint_file_name); // This is an initial condition, not a restart. Therefore, the current // time is set to 0, the iteration number as well and the pvd is cleared // from previous results - this->simulation_control->set_current_time(0.0); - this->simulation_control->set_iteration_number(0); - this->simulation_control->set_current_time_step( - this->simulation_parameters.simulation_control.dt); + // this->simulation_control->set_current_time(0.0); + // this->simulation_control->set_iteration_number(0); + // this->simulation_control->set_current_time_step( + // this->simulation_parameters.simulation_control.dt); - this->pvdhandler.times_and_names.clear(); + // this->pvdhandler.times_and_names.clear(); } else { diff --git a/source/solvers/navier_stokes_base.cc b/source/solvers/navier_stokes_base.cc index 073de2470e..8cb27fff42 100644 --- a/source/solvers/navier_stokes_base.cc +++ b/source/solvers/navier_stokes_base.cc @@ -1736,35 +1736,121 @@ void NavierStokesBase::read_checkpoint() { TimerOutput::Scope timer(this->computing_timer, "Read checkpoint"); + std::string prefix = this->simulation_parameters.simulation_control.output_folder + this->simulation_parameters.restart_parameters.filename; + + this->simulation_control->read(prefix); + this->pvdhandler.read(prefix); // When the average velocity profile is used as an initial condition, the // output folder path is not defined in the simulation control field but in // the initial condition field. - std::string checkpoint_folder; - std::string checkpoint_file_name; - if (this->simulation_parameters.initial_condition->type == - Parameters::InitialConditionType::average_velocity_profile && - !(this->simulation_parameters.restart_parameters.restart)) - { - checkpoint_folder = - this->simulation_parameters.initial_condition->checkpoint_folder; - checkpoint_file_name = - this->simulation_parameters.initial_condition->checkpoint_file_name; - } - else + // std::string checkpoint_folder; + // std::string checkpoint_file_name; + // if (this->simulation_parameters.initial_condition->type == + // Parameters::InitialConditionType::average_velocity_profile && + // !(this->simulation_parameters.restart_parameters.restart)) + // { + // checkpoint_folder = + // this->simulation_parameters.initial_condition->checkpoint_folder; + // checkpoint_file_name = + // this->simulation_parameters.initial_condition->checkpoint_file_name; + // } + // else + // { + // checkpoint_folder = + // this->simulation_parameters.simulation_control.output_folder; + // checkpoint_file_name = + // this->simulation_parameters.restart_parameters.filename; + // } + + this->pcout << "Reading checkpoint file: " << prefix << std::endl; + + this-> set_solution_from_checkpoint(prefix); + + // std::string prefix = checkpoint_folder + checkpoint_file_name; + + + if (simulation_parameters.flow_control.enable_flow_control) { - checkpoint_folder = - this->simulation_parameters.simulation_control.output_folder; - checkpoint_file_name = - this->simulation_parameters.restart_parameters.filename; + this->flow_control.read(prefix); } - std::string prefix = checkpoint_folder + checkpoint_file_name; + multiphysics->read_checkpoint(); - this->simulation_control->read(prefix); - this->pvdhandler.read(prefix); - const std::string filename = prefix + ".triangulation"; + // Deserialize all post-processing tables that are currently used + { + const Parameters::PostProcessing post_processing = + this->simulation_parameters.post_processing; + std::string prefix = + this->simulation_parameters.simulation_control.output_folder; + std::string suffix = ".checkpoint"; + if (post_processing.calculate_enstrophy) + deserialize_table(this->enstrophy_table, + prefix + post_processing.enstrophy_output_name + + suffix); + if (post_processing.calculate_pressure_power) + deserialize_table(this->pressure_power_table, + prefix + post_processing.pressure_power_output_name + + suffix); + if (post_processing.calculate_viscous_dissipation) + deserialize_table(this->viscous_dissipation_table, + prefix + + post_processing.viscous_dissipation_output_name + + suffix); + if (post_processing.calculate_kinetic_energy) + deserialize_table(this->kinetic_energy_table, + prefix + post_processing.kinetic_energy_output_name + + suffix); + if (post_processing.calculate_apparent_viscosity) + deserialize_table(this->apparent_viscosity_table, + prefix + + post_processing.apparent_viscosity_output_name + + suffix); + if (post_processing.calculate_flow_rate) + deserialize_table(this->flow_rate_table, + prefix + post_processing.flow_rate_output_name + + suffix); + if (post_processing.calculate_pressure_drop) + deserialize_table(this->pressure_drop_table, + prefix + post_processing.pressure_drop_output_name + + suffix); + if (this->simulation_parameters.forces_parameters.calculate_force) + for (unsigned int boundary_id = 0; + boundary_id < this->simulation_parameters.boundary_conditions.size; + ++boundary_id) + { + deserialize_table( + this->forces_tables[boundary_id], + prefix + + this->simulation_parameters.forces_parameters.force_output_name + + "_" + Utilities::int_to_string(boundary_id, 2) + suffix); + } + if (this->simulation_parameters.forces_parameters.calculate_torque) + for (unsigned int boundary_id = 0; + boundary_id < this->simulation_parameters.boundary_conditions.size; + ++boundary_id) + { + deserialize_table( + this->torques_tables[boundary_id], + prefix + + this->simulation_parameters.forces_parameters.torque_output_name + + "_" + Utilities::int_to_string(boundary_id, 2) + suffix); + } + if (this->simulation_parameters.analytical_solution->calculate_error()) + deserialize_table( + this->error_table, + prefix + + this->simulation_parameters.analytical_solution->get_filename() + + "_FD" + suffix); + } +} + +template +void +NavierStokesBase::set_solution_from_checkpoint(std::string checkpoint_file_prefix) +{ + const std::string filename = checkpoint_file_prefix + ".triangulation"; std::ifstream in(filename.c_str()); if (!in) AssertThrow(false, @@ -1777,8 +1863,10 @@ NavierStokesBase::read_checkpoint() try { if (auto tria = dynamic_cast *>( - this->triangulation.get())) - tria->load(filename.c_str()); + this->triangulation.get())){ + std::cout << filename << std::endl; + tria->load(filename.c_str()); + } } catch (...) { @@ -1786,6 +1874,8 @@ NavierStokesBase::read_checkpoint() ExcMessage("Cannot open snapshot mesh file or read the " "triangulation stored there.")); } + + setup_dofs(); enable_dynamic_zero_constraints_fd(); @@ -1817,12 +1907,10 @@ NavierStokesBase::read_checkpoint() parallel::distributed::SolutionTransfer system_trans_vectors( this->dof_handler); - if (simulation_parameters.post_processing.calculate_average_velocities || - this->simulation_parameters.initial_condition->type == - Parameters::InitialConditionType::average_velocity_profile) + if (simulation_parameters.post_processing.calculate_average_velocities || this->simulation_parameters.initial_condition->type == Parameters::InitialConditionType::average_velocity_profile) { std::vector sum_vectors = - this->average_velocities->read(prefix); + this->average_velocities->read(checkpoint_file_prefix); x_system.insert(x_system.end(), sum_vectors.begin(), sum_vectors.end()); } @@ -1834,9 +1922,7 @@ NavierStokesBase::read_checkpoint() previous_solutions[i] = distributed_previous_solutions[i]; } - if (simulation_parameters.post_processing.calculate_average_velocities || - this->simulation_parameters.initial_condition->type == - Parameters::InitialConditionType::average_velocity_profile) + if (simulation_parameters.post_processing.calculate_average_velocities || this->simulation_parameters.initial_condition->type == Parameters::InitialConditionType::average_velocity_profile) { this->average_velocities->calculate_average_velocities( this->local_evaluation_point, @@ -1845,94 +1931,6 @@ NavierStokesBase::read_checkpoint() simulation_control->get_time_step()); } - if (simulation_parameters.flow_control.enable_flow_control) - { - this->flow_control.read(prefix); - } - - // The average velocity profile initial condition uses the checkpoint - // mechanism to fetch a previous simulation solution. If an auxiliary physics - // was set to false in this previous simulation and is now set as true in a - // new simulation, the code would try to get an inexistant solution. - // Therefore, the multiphysics initial condition should not be set from the - // checkpoint. - if (this->simulation_parameters.initial_condition->type != - Parameters::InitialConditionType::average_velocity_profile || - (this->simulation_parameters.initial_condition->type == - Parameters::InitialConditionType::average_velocity_profile && - this->simulation_parameters.restart_parameters.restart)) - { - // std::cout << "Entering normally" << std::endl; - - multiphysics->read_checkpoint(); - } - - // Deserialize all post-processing tables that are currently used - { - const Parameters::PostProcessing post_processing = - this->simulation_parameters.post_processing; - std::string prefix = - this->simulation_parameters.simulation_control.output_folder; - std::string suffix = ".checkpoint"; - if (post_processing.calculate_enstrophy) - deserialize_table(this->enstrophy_table, - prefix + post_processing.enstrophy_output_name + - suffix); - if (post_processing.calculate_pressure_power) - deserialize_table(this->pressure_power_table, - prefix + post_processing.pressure_power_output_name + - suffix); - if (post_processing.calculate_viscous_dissipation) - deserialize_table(this->viscous_dissipation_table, - prefix + - post_processing.viscous_dissipation_output_name + - suffix); - if (post_processing.calculate_kinetic_energy) - deserialize_table(this->kinetic_energy_table, - prefix + post_processing.kinetic_energy_output_name + - suffix); - if (post_processing.calculate_apparent_viscosity) - deserialize_table(this->apparent_viscosity_table, - prefix + - post_processing.apparent_viscosity_output_name + - suffix); - if (post_processing.calculate_flow_rate) - deserialize_table(this->flow_rate_table, - prefix + post_processing.flow_rate_output_name + - suffix); - if (post_processing.calculate_pressure_drop) - deserialize_table(this->pressure_drop_table, - prefix + post_processing.pressure_drop_output_name + - suffix); - if (this->simulation_parameters.forces_parameters.calculate_force) - for (unsigned int boundary_id = 0; - boundary_id < this->simulation_parameters.boundary_conditions.size; - ++boundary_id) - { - deserialize_table( - this->forces_tables[boundary_id], - prefix + - this->simulation_parameters.forces_parameters.force_output_name + - "_" + Utilities::int_to_string(boundary_id, 2) + suffix); - } - if (this->simulation_parameters.forces_parameters.calculate_torque) - for (unsigned int boundary_id = 0; - boundary_id < this->simulation_parameters.boundary_conditions.size; - ++boundary_id) - { - deserialize_table( - this->torques_tables[boundary_id], - prefix + - this->simulation_parameters.forces_parameters.torque_output_name + - "_" + Utilities::int_to_string(boundary_id, 2) + suffix); - } - if (this->simulation_parameters.analytical_solution->calculate_error()) - deserialize_table( - this->error_table, - prefix + - this->simulation_parameters.analytical_solution->get_filename() + - "_FD" + suffix); - } } template