Example 9 - Solving a Transient Linear System in Parallel

This example shows how a simple, linear transient system can be solved in parallel. The system is simple scalar convection-diffusion with a specified external velocity. The initial condition is given, and the solution is advanced in time with a standard Crank-Nicholson time-stepping strategy.

C++ include files that we need

        #include <iostream>
        #include <algorithm>
        #include <math.h>

Basic include file needed for the mesh functionality.

        #include "libmesh.h"
        #include "mesh.h"
        #include "mesh_refinement.h"
        #include "gmv_io.h"
        #include "equation_systems.h"
        #include "fe.h"
        #include "quadrature_gauss.h"
        #include "dof_map.h"
        #include "sparse_matrix.h"
        #include "numeric_vector.h"
        #include "dense_matrix.h"
        #include "dense_vector.h"

Some (older) compilers do not offer full stream functionality, OStringStream works around this.

        #include "o_string_stream.h"

This example will solve a linear transient system, so we need to include the \p TransientLinearImplicitSystem definition.

        #include "linear_implicit_system.h"
        #include "transient_system.h"
        #include "vector_value.h"

The definition of a geometric element

        #include "elem.h"

Function prototype. This function will assemble the system matrix and right-hand-side at each time step. Note that since the system is linear we technically do not need to assmeble the matrix at each time step, but we will anyway. In subsequent examples we will employ adaptive mesh refinement, and with a changing mesh it will be necessary to rebuild the system matrix.

        void assemble_cd (EquationSystems& es,
        		  const std::string& system_name);

Function prototype. This function will initialize the system. Initialization functions are optional for systems. They allow you to specify the initial values of the solution. If an initialization function is not provided then the default (0) solution is provided.

        void init_cd (EquationSystems& es,
        	      const std::string& system_name);

Exact solution function prototype. This gives the exact solution as a function of space and time. In this case the initial condition will be taken as the exact solution at time 0, as will the Dirichlet boundary conditions at time t.

        Real exact_solution (const Real x,
        		     const Real y,
        		     const Real t);
        
        Number exact_value (const Point& p,
        		    const Parameters& parameters,
        		    const std::string&,
        		    const std::string&)
        {
          return exact_solution(p(0), p(1), parameters.get<Real> ("time"));
        }

We can now begin the main program. Note that this example will fail if you are using complex numbers since it was designed to be run only with real numbers.

        int main (int argc, char** argv)
        {

Initialize libMesh.

          libMesh::init (argc, argv);
        
        #ifndef ENABLE_AMR
          std::cerr << "ERROR: This example requires libMesh to be\n"
                    << "compiled with AMR support!"
                    << std::endl;
          return 0;
        #else
        
          {

Create a two-dimensional mesh.

            Mesh mesh (2);

Read the mesh from file. This is the coarse mesh that will be used in example 10 to demonstrate adaptive mesh refinement. Here we will simply read it in and uniformly refine it 5 times before we compute with it.

            mesh.read ("../ex10/mesh.xda");

Create a MeshRefinement object to handle refinement of our mesh. This class handles all the details of mesh refinement and coarsening.

            MeshRefinement mesh_refinement (mesh);

Uniformly refine the mesh 5 times. This is the first time we use the mesh refinement capabilities of the library.

            mesh_refinement.uniformly_refine (5);

Print information about the mesh to the screen.

            mesh.print_info();

Create an equation systems object.

            EquationSystems equation_systems (mesh);

Add a transient system to the EquationSystems object named "Convection-Diffusion".

            TransientLinearImplicitSystem & system = 
              equation_systems.add_system<TransientLinearImplicitSystem> ("Convection-Diffusion");

Adds the variable "u" to "Convection-Diffusion". "u" will be approximated using first-order approximation.

            system.add_variable ("u", FIRST);

Give the system a pointer to the matrix assembly and initialization functions.

            system.attach_assemble_function (assemble_cd);
            system.attach_init_function (init_cd);

Initialize the data structures for the equation system.

            equation_systems.init ();

Prints information about the system to the screen.

            equation_systems.print_info();

Write out the initial conditions.

            GMVIO(mesh).write_equation_systems ("out_000.gmv",
        					equation_systems);

The Convection-Diffusion system requires that we specify the flow velocity. We will specify it as a RealVectorValue data type and then use the Parameters object to pass it to the assemble function.

            equation_systems.parameters.set<RealVectorValue>("velocity") = 
              RealVectorValue (0.8, 0.8);

Solve the system "Convection-Diffusion". This will be done by looping over the specified time interval and calling the solve() member at each time step. This will assemble the system and call the linear solver.

            const Real dt = 0.025;
            Real time     = 0.;
            
            for (unsigned int t_step = 0; t_step < 50; t_step++)
              {

Incremenet the time counter, set the time and the time step size as parameters in the EquationSystem.

                time += dt;
        
        	equation_systems.parameters.set<Real> ("time") = time;
        	equation_systems.parameters.set<Real> ("dt")   = dt;

A pretty update message

                std::cout << " Solving time step ";

Since some compilers fail to offer full stream functionality, libMesh offers a string stream to work around this. Note that for other compilers, this is just a set of preprocessor macros and therefore should cost nothing (compared to a hand-coded string stream). We use additional curly braces here simply to enforce data locality.

                {
        	  OStringStream out;
        
        	  OSSInt(out,2,t_step);
        	  out << ", time=";
        	  OSSRealzeroleft(out,6,3,time);
        	  out <<  "...";
        	  std::cout << out.str() << std::endl;
        	}

At this point we need to update the old solution vector. The old solution vector will be the current solution vector from the previous time step. We will do this by extracting the system from the \p EquationSystems object and using vector assignment. Since only \p TransientSystems (and systems derived from them) contain old solutions we need to specify the system type when we ask for it.

                TransientLinearImplicitSystem&  system =
        	  equation_systems.get_system<TransientLinearImplicitSystem>("Convection-Diffusion");
        
        	*system.old_local_solution = *system.current_local_solution;

Assemble & solve the linear system

                equation_systems.get_system("Convection-Diffusion").solve();

Output evey 10 timesteps to file.

                if ( (t_step+1)%10 == 0)
        	  {
        	    OStringStream file_name;
        
        	    file_name << "out_";
        	    OSSRealzeroright(file_name,3,0,t_step+1);
        	    file_name << ".gmv";
        
        	    GMVIO(mesh).write_equation_systems (file_name.str(),
        						equation_systems);
        	  }
              }
          }
        #endif // #ifdef ENABLE_AMR

All done.

          return libMesh::close ();
        }

We now define the function which provides the initialization routines for the "Convection-Diffusion" system. This handles things like setting initial conditions and boundary conditions.

        void init_cd (EquationSystems& es,
        	      const std::string& system_name)
        {

It is a good idea to make sure we are initializing the proper system.

          assert (system_name == "Convection-Diffusion");

Get a reference to the Convection-Diffusion system object.

          TransientLinearImplicitSystem & system =
            es.get_system<TransientLinearImplicitSystem>("Convection-Diffusion");

Project initial conditions at time 0

          es.parameters.set<Real> ("time") = 0;
          
          system.project_solution(exact_value, NULL, es.parameters);
        }

Now we define the assemble function which will be used by the EquationSystems object at each timestep to assemble the linear system for solution.

        void assemble_cd (EquationSystems& es,
        		  const std::string& system_name)
        {
        #ifdef ENABLE_AMR

It is a good idea to make sure we are assembling the proper system.

          assert (system_name == "Convection-Diffusion");

Get a constant reference to the mesh object.

          const Mesh& mesh = es.get_mesh();

The dimension that we are running

          const unsigned int dim = mesh.mesh_dimension();

Get a reference to the Convection-Diffusion system object.

          TransientLinearImplicitSystem & system =
            es.get_system<TransientLinearImplicitSystem> ("Convection-Diffusion");

Get a constant reference to the Finite Element type for the first (and only) variable in the system.

          FEType fe_type = system.variable_type(0);

Build a Finite Element object of the specified type. Since the \p FEBase::build() member dynamically creates memory we will store the object as an \p AutoPtr. This can be thought of as a pointer that will clean up after itself.

          AutoPtr<FEBase> fe      (FEBase::build(dim, fe_type));
          AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));

A Gauss quadrature rule for numerical integration. Let the \p FEType object decide what order rule is appropriate.

          QGauss qrule (dim,   fe_type.default_quadrature_order());
          QGauss qface (dim-1, fe_type.default_quadrature_order());

Tell the finite element object to use our quadrature rule.

          fe->attach_quadrature_rule      (&qrule);
          fe_face->attach_quadrature_rule (&qface);

Here we define some references to cell-specific data that will be used to assemble the linear system. We will start with the element Jacobian * quadrature weight at each integration point.

          const std::vector<Real>& JxW      = fe->get_JxW();
          const std::vector<Real>& JxW_face = fe_face->get_JxW();

The element shape functions evaluated at the quadrature points.

          const std::vector<std::vector<Real> >& phi = fe->get_phi();
          const std::vector<std::vector<Real> >& psi = fe_face->get_phi();

The element shape function gradients evaluated at the quadrature points.

          const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();

The XY locations of the quadrature points used for face integration

          const std::vector<Point>& qface_points = fe_face->get_xyz();

A reference to the \p DofMap object for this system. The \p DofMap object handles the index translation from node and element numbers to degree of freedom numbers. We will talk more about the \p DofMap in future examples.

          const DofMap& dof_map = system.get_dof_map();

Define data structures to contain the element matrix and right-hand-side vector contribution. Following basic finite element terminology we will denote these "Ke" and "Fe".

          DenseMatrix<Number> Ke;
          DenseVector<Number> Fe;

This vector will hold the degree of freedom indices for the element. These define where in the global system the element degrees of freedom get mapped.

          std::vector<unsigned int> dof_indices;

Here we extract the velocity & parameters that we put in the EquationSystems object.

          const RealVectorValue velocity =
            es.parameters.get<RealVectorValue> ("velocity");
        
          const Real dt = es.parameters.get<Real>   ("dt");
          const Real time = es.parameters.get<Real> ("time");

Now we will loop over all the elements in the mesh that live on the local processor. We will compute the element matrix and right-hand-side contribution. Since the mesh will be refined we want to only consider the ACTIVE elements, hence we use a variant of the \p active_elem_iterator. const_active_local_elem_iterator el (mesh.elements_begin()); const const_active_local_elem_iterator end_el (mesh.elements_end());

          MeshBase::const_element_iterator       el     = mesh.active_local_elements_begin();
          const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end(); 
        
          for ( ; el != end_el; ++el)
            {

Store a pointer to the element we are currently working on. This allows for nicer syntax later.

              const Elem* elem = *el;

Get the degree of freedom indices for the current element. These define where in the global matrix and right-hand-side this element will contribute to.

              dof_map.dof_indices (elem, dof_indices);

Compute the element-specific data for the current element. This involves computing the location of the quadrature points (q_point) and the shape functions (phi, dphi) for the current element.

              fe->reinit (elem);

Zero the element matrix and right-hand side before summing them. We use the resize member here because the number of degrees of freedom might have changed from the last element. Note that this will be the case if the element type is different (i.e. the last element was a triangle, now we are on a quadrilateral).

              Ke.resize (dof_indices.size(),
        		 dof_indices.size());
        
              Fe.resize (dof_indices.size());

Now we will build the element matrix and right-hand-side. Constructing the RHS requires the solution and its gradient from the previous timestep. This myst be calculated at each quadrature point by summing the solution degree-of-freedom values by the appropriate weight functions.

              for (unsigned int qp=0; qp<qrule.n_points(); qp++)
        	{

Values to hold the old solution & its gradient.

                  Number   u_old = 0.;
        	  Gradient grad_u_old;

Compute the old solution & its gradient.

                  for (unsigned int l=0; l<phi.size(); l++)
        	    {
        	      u_old      += phi[l][qp]*system.old_solution  (dof_indices[l]);

This will work, grad_u_old += dphi[l][qp]*system.old_solution (dof_indices[l]); but we can do it without creating a temporary like this:

                      grad_u_old.add_scaled (dphi[l][qp],system.old_solution (dof_indices[l]));
        	    }

Now compute the element matrix and RHS contributions.

                  for (unsigned int i=0; i<phi.size(); i++)
        	    {

The RHS contribution

                      Fe(i) += JxW[qp]*(

Mass matrix term

                                        u_old*phi[i][qp] +
        				-.5*dt*(

Convection term (grad_u_old may be complex, so the order here is important!)

                                                (grad_u_old*velocity)*phi[i][qp] +

Diffusion term

                                                0.01*(grad_u_old*dphi[i][qp]))     
        				);
        	      
        	      for (unsigned int j=0; j<phi.size(); j++)
        		{

The matrix contribution

                          Ke(i,j) += JxW[qp]*(

Mass-matrix

                                              phi[i][qp]*phi[j][qp] + 
        
        				      .5*dt*(

Convection term

                                                     (velocity*dphi[j][qp])*phi[i][qp] +

Diffusion term

                                                     0.01*(dphi[i][qp]*dphi[j][qp]))      
        				      );
        		} 
        	    } 
        	}

At this point the interior element integration has been completed. However, we have not yet addressed boundary conditions. For this example we will only consider simple Dirichlet boundary conditions imposed via the penalty method.

The following loops over the sides of the element. If the element has no neighbor on a side then that side MUST live on a boundary of the domain.

The penalty value.

                const Real penalty = 1.e10;

The following loops over the sides of the element. If the element has no neighbor on a side then that side MUST live on a boundary of the domain.

                for (unsigned int s=0; s<elem->n_sides(); s++)
        	  if (elem->neighbor(s) == NULL)
        	    {
        	      fe_face->reinit(elem,s);
        	      
        	      for (unsigned int qp=0; qp<qface.n_points(); qp++)
        		{
        		  const Number value = exact_solution (qface_points[qp](0),
        						       qface_points[qp](1),
        						       time);

RHS contribution

                          for (unsigned int i=0; i<psi.size(); i++)
        		    Fe(i) += penalty*JxW_face[qp]*value*psi[i][qp];

Matrix contribution

                          for (unsigned int i=0; i<psi.size(); i++)
        		    for (unsigned int j=0; j<psi.size(); j++)
        		      Ke(i,j) += penalty*JxW_face[qp]*psi[i][qp]*psi[j][qp];
        		}
        	    } 
              }

The element matrix and right-hand-side are now built for this element. Add them to the global matrix and right-hand-side vector. The \p PetscMatrix::add_matrix() and \p PetscVector::add_vector() members do this for us.

              system.matrix->add_matrix (Ke, dof_indices);
              system.rhs->add_vector    (Fe, dof_indices);
            }

That concludes the system matrix assembly routine.

        #endif // #ifdef ENABLE_AMR
        }

The program without comments:

 
  
  #include <iostream>
  #include <algorithm>
  #include <math.h>
  
  #include "libmesh.h"
  #include "mesh.h"
  #include "mesh_refinement.h"
  #include "gmv_io.h"
  #include "equation_systems.h"
  #include "fe.h"
  #include "quadrature_gauss.h"
  #include "dof_map.h"
  #include "sparse_matrix.h"
  #include "numeric_vector.h"
  #include "dense_matrix.h"
  #include "dense_vector.h"
  
  #include "o_string_stream.h"
  
  #include "linear_implicit_system.h"
  #include "transient_system.h"
  #include "vector_value.h"
  
  #include "elem.h"
  
  void assemble_cd (EquationSystems& es,
  		  const std::string& system_name);
  
  void init_cd (EquationSystems& es,
  	      const std::string& system_name);
  
  Real exact_solution (const Real x,
  		     const Real y,
  		     const Real t);
  
  Number exact_value (const Point& p,
  		    const Parameters& parameters,
  		    const std::string&,
  		    const std::string&)
  {
    return exact_solution(p(0), p(1), parameters.get<Real> ("time"));
  }
  
  
  
  int main (int argc, char** argv)
  {
    libMesh::init (argc, argv);
  
  #ifndef ENABLE_AMR
    std::cerr << "ERROR: This example requires libMesh to be\n"
              << "compiled with AMR support!"
              << std::endl;
    return 0;
  #else
  
    {    
      Mesh mesh (2);
          
      mesh.read ("../ex10/mesh.xda");
      
      MeshRefinement mesh_refinement (mesh);
      
      mesh_refinement.uniformly_refine (5);
      
      mesh.print_info();
      
      EquationSystems equation_systems (mesh);
      
      TransientLinearImplicitSystem & system = 
        equation_systems.add_system<TransientLinearImplicitSystem> ("Convection-Diffusion");
        
      system.add_variable ("u", FIRST);
        
      system.attach_assemble_function (assemble_cd);
      system.attach_init_function (init_cd);
        
      equation_systems.init ();
        
      equation_systems.print_info();
        
      GMVIO(mesh).write_equation_systems ("out_000.gmv",
  					equation_systems);
      
      equation_systems.parameters.set<RealVectorValue>("velocity") = 
        RealVectorValue (0.8, 0.8);
      
      const Real dt = 0.025;
      Real time     = 0.;
      
      for (unsigned int t_step = 0; t_step < 50; t_step++)
        {
  	time += dt;
  
  	equation_systems.parameters.set<Real> ("time") = time;
  	equation_systems.parameters.set<Real> ("dt")   = dt;
  
  	std::cout << " Solving time step ";
  	
  	{
  	  OStringStream out;
  
  	  OSSInt(out,2,t_step);
  	  out << ", time=";
  	  OSSRealzeroleft(out,6,3,time);
  	  out <<  "...";
  	  std::cout << out.str() << std::endl;
  	}
  	
  	TransientLinearImplicitSystem&  system =
  	  equation_systems.get_system<TransientLinearImplicitSystem>("Convection-Diffusion");
  
  	*system.old_local_solution = *system.current_local_solution;
  	
  	equation_systems.get_system("Convection-Diffusion").solve();
  	
  	if ( (t_step+1)%10 == 0)
  	  {
  	    OStringStream file_name;
  
  	    file_name << "out_";
  	    OSSRealzeroright(file_name,3,0,t_step+1);
  	    file_name << ".gmv";
  
  	    GMVIO(mesh).write_equation_systems (file_name.str(),
  						equation_systems);
  	  }
        }
    }
  #endif // #ifdef ENABLE_AMR
  
    return libMesh::close ();
  }
  
  void init_cd (EquationSystems& es,
  	      const std::string& system_name)
  {
    assert (system_name == "Convection-Diffusion");
  
    TransientLinearImplicitSystem & system =
      es.get_system<TransientLinearImplicitSystem>("Convection-Diffusion");
  
    es.parameters.set<Real> ("time") = 0;
    
    system.project_solution(exact_value, NULL, es.parameters);
  }
  
  
  
  void assemble_cd (EquationSystems& es,
  		  const std::string& system_name)
  {
  #ifdef ENABLE_AMR
    assert (system_name == "Convection-Diffusion");
    
    const Mesh& mesh = es.get_mesh();
    
    const unsigned int dim = mesh.mesh_dimension();
    
    TransientLinearImplicitSystem & system =
      es.get_system<TransientLinearImplicitSystem> ("Convection-Diffusion");
    
    FEType fe_type = system.variable_type(0);
    
    AutoPtr<FEBase> fe      (FEBase::build(dim, fe_type));
    AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
    
    QGauss qrule (dim,   fe_type.default_quadrature_order());
    QGauss qface (dim-1, fe_type.default_quadrature_order());
  
    fe->attach_quadrature_rule      (&qrule);
    fe_face->attach_quadrature_rule (&qface);
  
    const std::vector<Real>& JxW      = fe->get_JxW();
    const std::vector<Real>& JxW_face = fe_face->get_JxW();
    
    const std::vector<std::vector<Real> >& phi = fe->get_phi();
    const std::vector<std::vector<Real> >& psi = fe_face->get_phi();
  
    const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
  
    const std::vector<Point>& qface_points = fe_face->get_xyz();
    
    const DofMap& dof_map = system.get_dof_map();
  
    DenseMatrix<Number> Ke;
    DenseVector<Number> Fe;
    
    std::vector<unsigned int> dof_indices;
  
    const RealVectorValue velocity =
      es.parameters.get<RealVectorValue> ("velocity");
  
    const Real dt = es.parameters.get<Real>   ("dt");
    const Real time = es.parameters.get<Real> ("time");
  
  
    MeshBase::const_element_iterator       el     = mesh.active_local_elements_begin();
    const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end(); 
  
    for ( ; el != end_el; ++el)
      {    
        const Elem* elem = *el;
        
        dof_map.dof_indices (elem, dof_indices);
        
        fe->reinit (elem);
        
        Ke.resize (dof_indices.size(),
  		 dof_indices.size());
  
        Fe.resize (dof_indices.size());
        
        for (unsigned int qp=0; qp<qrule.n_points(); qp++)
  	{
  	  Number   u_old = 0.;
  	  Gradient grad_u_old;
  	  
  	  for (unsigned int l=0; l<phi.size(); l++)
  	    {
  	      u_old      += phi[l][qp]*system.old_solution  (dof_indices[l]);
  	      
  	      grad_u_old.add_scaled (dphi[l][qp],system.old_solution (dof_indices[l]));
  	    }
  	  
  	  for (unsigned int i=0; i<phi.size(); i++)
  	    {
  	      Fe(i) += JxW[qp]*(
  				u_old*phi[i][qp] +
  				-.5*dt*(
  					(grad_u_old*velocity)*phi[i][qp] +
  					
  					0.01*(grad_u_old*dphi[i][qp]))     
  				);
  	      
  	      for (unsigned int j=0; j<phi.size(); j++)
  		{
  		  Ke(i,j) += JxW[qp]*(
  				      phi[i][qp]*phi[j][qp] + 
  
  				      .5*dt*(
  					     (velocity*dphi[j][qp])*phi[i][qp] +
  					     
  					     0.01*(dphi[i][qp]*dphi[j][qp]))      
  				      );
  		} 
  	    } 
  	} 
  
        {
  	const Real penalty = 1.e10;
  
  	for (unsigned int s=0; s<elem->n_sides(); s++)
  	  if (elem->neighbor(s) == NULL)
  	    {
  	      fe_face->reinit(elem,s);
  	      
  	      for (unsigned int qp=0; qp<qface.n_points(); qp++)
  		{
  		  const Number value = exact_solution (qface_points[qp](0),
  						       qface_points[qp](1),
  						       time);
  						       
  		  for (unsigned int i=0; i<psi.size(); i++)
  		    Fe(i) += penalty*JxW_face[qp]*value*psi[i][qp];
  
  		  for (unsigned int i=0; i<psi.size(); i++)
  		    for (unsigned int j=0; j<psi.size(); j++)
  		      Ke(i,j) += penalty*JxW_face[qp]*psi[i][qp]*psi[j][qp];
  		}
  	    } 
        } 
  
        system.matrix->add_matrix (Ke, dof_indices);
        system.rhs->add_vector    (Fe, dof_indices);
      }
    
  #endif // #ifdef ENABLE_AMR
  }

The console output of the program:

***************************************************************
* Running Example  ./ex9-devel
***************************************************************
 
 Mesh Information:
  mesh_dimension()=2
  spatial_dimension()=3
  n_nodes()=6273
  n_elem()=13650
   n_local_elem()=13650
   n_active_elem()=10240
  n_subdomains()=1
  n_processors()=1
  processor_id()=0

 EquationSystems
  n_systems()=1
   System "Convection-Diffusion"
    Type "TransientLinearImplicit"
    Variables="u" 
    Finite Element Types="LAGRANGE" 
    Approximation Orders="FIRST" 
    n_dofs()=6273
    n_local_dofs()=6273
    n_constrained_dofs()=0
    n_vectors()=3

 Solving time step  0, time=0.0250...
 Solving time step  1, time=0.0500...
 Solving time step  2, time=0.0750...
 Solving time step  3, time=0.1000...
 Solving time step  4, time=0.1250...
 Solving time step  5, time=0.1500...
 Solving time step  6, time=0.1750...
 Solving time step  7, time=0.2000...
 Solving time step  8, time=0.2250...
 Solving time step  9, time=0.2500...
 Solving time step 10, time=0.2750...
 Solving time step 11, time=0.3000...
 Solving time step 12, time=0.3250...
 Solving time step 13, time=0.3500...
 Solving time step 14, time=0.3750...
 Solving time step 15, time=0.4000...
 Solving time step 16, time=0.4250...
 Solving time step 17, time=0.4500...
 Solving time step 18, time=0.4750...
 Solving time step 19, time=0.5000...
 Solving time step 20, time=0.5250...
 Solving time step 21, time=0.5500...
 Solving time step 22, time=0.5750...
 Solving time step 23, time=0.6000...
 Solving time step 24, time=0.6250...
 Solving time step 25, time=0.6500...
 Solving time step 26, time=0.6750...
 Solving time step 27, time=0.7000...
 Solving time step 28, time=0.7250...
 Solving time step 29, time=0.7500...
 Solving time step 30, time=0.7750...
 Solving time step 31, time=0.8000...
 Solving time step 32, time=0.8250...
 Solving time step 33, time=0.8500...
 Solving time step 34, time=0.8750...
 Solving time step 35, time=0.9000...
 Solving time step 36, time=0.9250...
 Solving time step 37, time=0.9500...
 Solving time step 38, time=0.9750...
 Solving time step 39, time=1.0000...
 Solving time step 40, time=1.0300...
 Solving time step 41, time=1.0500...
 Solving time step 42, time=1.0800...
 Solving time step 43, time=1.1000...
 Solving time step 44, time=1.1200...
 Solving time step 45, time=1.1500...
 Solving time step 46, time=1.1700...
 Solving time step 47, time=1.2000...
 Solving time step 48, time=1.2200...
 Solving time step 49, time=1.2500...
 
***************************************************************
* Done Running Example  ./ex9-devel
***************************************************************