Skip to content

Commit

Permalink
added set_solution_from_checkpoint function to clean the use of the r…
Browse files Browse the repository at this point in the history 
…estart with the initial condition. Scaling bug in average velocity
  • Loading branch information
@mivaia
mivaia committed Aug 2, 2024
1 parent 7ae0ec2 commit fbe5727
Show file tree
Hide file tree
Showing 3 changed files with 140 additions and 131 deletions.
8 changes: 7 additions & 1 deletion include/solvers/navier_stokes_base.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -382,11 +382,17 @@ class NavierStokesBase : public PhysicsSolver<VectorType>
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
Expand Down
31 changes: 18 additions & 13 deletions source/solvers/gls_navier_stokes.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1281,29 +1281,34 @@ GLSNavierStokesSolver<dim>::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
{
Expand Down
232 changes: 115 additions & 117 deletions source/solvers/navier_stokes_base.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1736,35 +1736,121 @@ void
NavierStokesBase<dim, VectorType, DofsType>::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 <int dim, typename VectorType, typename DofsType>
void
NavierStokesBase<dim, VectorType, DofsType>::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,
Expand All @@ -1777,15 +1863,19 @@ NavierStokesBase<dim, VectorType, DofsType>::read_checkpoint()
try
{
if (auto tria = dynamic_cast<parallel::distributed::Triangulation<dim> *>(
this->triangulation.get()))
tria->load(filename.c_str());
this->triangulation.get())){
std::cout << filename << std::endl;
tria->load(filename.c_str());
}
}
catch (...)
{
AssertThrow(false,
ExcMessage("Cannot open snapshot mesh file or read the "
"triangulation stored there."));
}


setup_dofs();
enable_dynamic_zero_constraints_fd();

Expand Down Expand Up @@ -1817,12 +1907,10 @@ NavierStokesBase<dim, VectorType, DofsType>::read_checkpoint()
parallel::distributed::SolutionTransfer<dim, VectorType> 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<VectorType *> 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());
}
Expand All @@ -1834,9 +1922,7 @@ NavierStokesBase<dim, VectorType, DofsType>::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,
Expand All @@ -1845,94 +1931,6 @@ NavierStokesBase<dim, VectorType, DofsType>::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 <int dim, typename VectorType, typename DofsType>
Expand Down

0 comments on commit fbe5727

Please sign in to comment.