oomph-lib

Simulations of time-dependent problem can be very time consuming and it is important to be able to restart simulations, e.g. to continue a run after a system crash, etc. We shall illustrate oomph-lib's dump/restart capabilities by re-visiting the 2D unsteady heat equation discussed in a previous example.

As before, we consider the unforced case, $ f=0 $

, and choose boundary and initial conditions that are consistent with the exact solution

$u_0(x_1,x_2,t) = e^{-Kt}\sin\left( \sqrt{K} \left( x_1 \cos \Phi + x_2 \sin \Phi\right)\right), \ \ \ \ \ \ \ \ \ \ (4)$

where

and $\Phi$ are constants, controlling the decay rate of the solution and its spatial orientation, respectively.

The figure below shows a plot of computed and exact solution at a control node as a function of time. The solid lines represent quantitites computed during the original simulation; the dashed line shows the corresponding data from a second simulation that was restarted with the restart file generated at timestep 23 of the original simulation.

Most of the driver code for this example is identical to that discussed in the previous example, therefore we only discuss the modifications required to enable the dump and restart operations:

Global parameters and functions

The driver code

//=======start_of_main====================================================
/// \short Driver code for unsteady heat equation with option for
/// restart from disk: Only a single command line argument is allowed.
/// If specified it is interpreted as the name of the restart file.
//========================================================================
int main(int argc, char* argv[])
{

 // Store command line arguments
 CommandLineArgs::setup(argc,argv);

The problem class

 /// \short Dump problem to disk to allow for restart.
 void dump_it(ofstream& dump_file);

 /// \short Read problem for restart from specified restart file.
 void restart(ifstream& restart_file);

The problem constructor

The problem destructor

Actions before timestep

Set initial condition

//======================start_of_set_initial_condition====================
/// \short Set initial condition: Assign previous and current values
/// from exact solution or from restart file.
//========================================================================
template<class ELEMENT>
void UnsteadyHeatProblem<ELEMENT>::set_initial_condition()
{ 

 // Pointer to restart file
 ifstream* restart_file_pt=0;

 // Restart?
 //---------
 // Restart file specified via command line [all programs have at least
 // a single command line argument: their name. Ignore this here.]
 if (CommandLineArgs::Argc==1)
  {
   cout << "No restart -- setting IC from exact solution" << std::endl;
  }
 else if (CommandLineArgs::Argc==2)
  {
   // Open restart file
   restart_file_pt= new ifstream(CommandLineArgs::Argv[1],ios_base::in);
   if (restart_file_pt!=0)
    {
     cout << "Have opened " << CommandLineArgs::Argv[1] << 
      " for restart. " << std::endl;
    }
   else
    {
     std::ostringstream error_stream;
     error_stream 
      << "ERROR while trying to open " << CommandLineArgs::Argv[1] << 
      " for restart." << std::endl;

     throw OomphLibError(
      error_stream.str(),
      "UnsteadyHeatProblem<ELEMENT>::set_initial_condition()",
      OOMPH_EXCEPTION_LOCATION);
    }
  }
 // More than one command line argument?
 else 
  {
   std::ostringstream error_stream;
   error_stream << "Can only specify one input file\n" 
                << "You specified the following command line arguments:\n";
   //Fix this
   CommandLineArgs::output();

   throw OomphLibError( 
    error_stream.str(),
    "UnsteadyHeatProblem<ELEMENT>::set_initial_condition()",
    OOMPH_EXCEPTION_LOCATION);
  }


 // Read restart data:
 //-------------------
 if (restart_file_pt!=0)
  {

   // Read the problem data from the restart file
   restart(*restart_file_pt);

  }
 // Assign initial condition from exact solution
 //---------------------------------------------
 else
  {
   
   // Choose timestep
   double dt=0.01;
   
   // Initialise timestep -- also sets the weights for all timesteppers
   // in the problem.
   initialise_dt(dt);
   
   // Backup time in global Time object
   double backed_up_time=time_pt()->time();
   
   // Past history fo needs to be established for t=time0-deltat, ...
   // Then provide current values (at t=time0) which will also form
   // the initial guess for the first solve at t=time0+deltat
   
   // Vector of exact solution value
   Vector<double> soln(1);
   Vector<double> x(2);
   
   //Find number of nodes in mesh
   unsigned num_nod = mesh_pt()->nnode();
   
   // Set continuous times at previous timesteps
   int nprev_steps=time_stepper_pt()->nprev_values();
   Vector<double> prev_time(nprev_steps+1);
   for (int itime=nprev_steps;itime>=0;itime--)
    {
     prev_time[itime]=time_pt()->time(unsigned(itime));
    } 
   
   // Loop over current & previous timesteps
   for (int itime=nprev_steps;itime>=0;itime--)
    {
     // Continuous time
     double time=prev_time[itime];
     cout << "setting IC at time =" << time << std::endl;
     
     // Loop over the nodes to set initial guess everywhere
     for (unsigned n=0;n<num_nod;n++)
      {
       // Get nodal coordinates
       x[0]=mesh_pt()->node_pt(n)->x(0);
       x[1]=mesh_pt()->node_pt(n)->x(1);
       
       // Get intial solution
       ExactSolnForUnsteadyHeat::get_exact_u(time,x,soln);
       
       // Assign solution
       mesh_pt()->node_pt(n)->set_value(itime,0,soln[0]);
       
       // Loop over coordinate directions: Previous position = present position
       for (unsigned i=0;i<2;i++)
        {
         mesh_pt()->node_pt(n)->x(itime,i)=x[i];
        }
      } 
    }
   
   // Reset backed up time for global timestepper
   time_pt()->time()=backed_up_time;

  }


} // end of set_initial_condition

Post-processing

 // Write restart file
 sprintf(filename,"%s/restart%i.dat",doc_info.directory().c_str(),
         doc_info.number());
 some_file.open(filename);
 dump_it(some_file);
 some_file.close();

Dumping the solution

In the present problem, the generic Problem::dump(...) function is sufficient to record the current state of the simulation, therefore no additional information needs to be added to the dump file. The section Comments and Exercises below contains an exercise that illustrates how to customise the dump function to record additional parameters; the demo code with spatial adaptivity provides another example of a customised dump/restart function.

Example problem: Solution of the 2D unsteady heat equation with restarts

Global parameters and functions

The driver code

The problem class

The problem constructor

The problem destructor

Actions before timestep

Set initial condition

Post-processing

Dumping the solution

Reading a solution from disk

Comments and Exercises

Exercises

Source files for this tutorial

PDF file