Example 13 - Unsteady Navier-Stokes Equations - Unsteady Nonlinear Systems of Equations

This example shows how a simple, unsteady, nonlinear system of equations can be solved in parallel. The system of equations are the familiar Navier-Stokes equations for low-speed incompressible fluid flow. This example introduces the concept of the inner nonlinear loop for each timestep, and requires a good deal of linear algebra number-crunching at each step. If you have the General Mesh Viewer (GMV) installed, the script movie.sh in this directory will also take appropriate screen shots of each of the solution files in the time sequence. These rgb files can then be animated with the "animate" utility of ImageMagick if it is installed on your system. On a PIII 1GHz machine in debug mode, this example takes a little over a minute to run. If you would like to see a more detailed time history, or compute more timesteps, that is certainly possible by changing the n_timesteps and dt variables below.

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_generation.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 "linear_implicit_system.h"
        #include "transient_system.h"
        #include "perf_log.h"
        #include "boundary_info.h"
        #include "utility.h"

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

        #include "o_string_stream.h"

For systems of equations the \p DenseSubMatrix and \p DenseSubVector provide convenient ways for assembling the element matrix and vector on a component-by-component basis.

        #include "dense_submatrix.h"
        #include "dense_subvector.h"

The definition of a geometric element

        #include "elem.h"

Function prototype. This function will assemble the system matrix and right-hand-side.

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

The main program.

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

Initialize libMesh.

          libMesh::init (argc, argv);
          {

Set the dimensionality of the mesh = 2

            const unsigned int dim = 2;

Create a two-dimensional mesh.

            Mesh mesh (dim);

Use the MeshTools::Generation mesh generator to create a uniform grid on the square [-1,1]^D. We instruct the mesh generator to build a mesh of 8x8 \p Quad9 elements in 2D, or \p Hex27 elements in 3D. Building these higher-order elements allows us to use higher-order approximation, as in example 3.

            MeshTools::Generation::build_square (mesh,
        					 20, 20,
        					 0., 1.,
        					 0., 1.,
        					 QUAD9);

Print information about the mesh to the screen.

            mesh.print_info();

Create an equation systems object.

            EquationSystems equation_systems (mesh);

Declare the system and its variables.

Creates a transient system named "Navier-Stokes"

              TransientLinearImplicitSystem & system = 
        	equation_systems.add_system<TransientLinearImplicitSystem> ("Navier-Stokes");

Add the variables "u" & "v" to "Navier-Stokes". They will be approximated using second-order approximation.

              system.add_variable ("u", SECOND);
              system.add_variable ("v", SECOND);

Add the variable "p" to "Navier-Stokes". This will be approximated with a first-order basis, providing an LBB-stable pressure-velocity pair.

              system.add_variable ("p", FIRST);

Give the system a pointer to the matrix assembly function.

              system.attach_assemble_function (assemble_stokes);

Initialize the data structures for the equation system.

              equation_systems.init ();

Prints information about the system to the screen.

              equation_systems.print_info();
            }

Create a performance-logging object for this example

            PerfLog perf_log("Example 13");

Now we begin the timestep loop to compute the time-accurate solution of the equations.

            const Real dt = 0.01;
            Real time     = 0.0;
            const unsigned int n_timesteps = 15;

The number of steps and the stopping criterion are also required for the nonlinear iterations.

            const unsigned int n_nonlinear_steps = 15;
            const Real nonlinear_tolerance       = 1.e-3;

We also set a standard linear solver flag in the EquationSystems object which controls the maxiumum number of linear solver iterations allowed.

            equation_systems.parameters.set<unsigned int>("linear solver maximum iterations") = 250;

Tell the system of equations what the timestep is by using the set_parameter function. The matrix assembly routine can then reference this parameter.

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

Get a reference to the Stokes system to use later.

            TransientLinearImplicitSystem&  navier_stokes_system =
        	  equation_systems.get_system<TransientLinearImplicitSystem>("Navier-Stokes");

The first thing to do is to get a copy of the solution at the current nonlinear iteration. This value will be used to determine if we can exit the nonlinear loop.

            AutoPtr<NumericVector<Number> >
              last_nonlinear_soln (navier_stokes_system.solution->clone());
        
            for (unsigned int t_step=0; t_step<n_timesteps; ++t_step)
              {

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

                time += dt;

Let the system of equations know the current time. This might be necessary for a time-dependent forcing function for example.

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

A pretty update message

                std::cout << "\n\n*** Solving time step " << t_step << ", time = " << time << " ***" << std::endl;

Now we need to update the solution vector from the previous time step. This is done directly through the reference to the Stokes system.

                *navier_stokes_system.old_local_solution = *navier_stokes_system.current_local_solution;

At the beginning of each solve, reset the linear solver tolerance to a "reasonable" starting value.

                const Real initial_linear_solver_tol = 1.e-6;
        	equation_systems.parameters.set<Real> ("linear solver tolerance") = initial_linear_solver_tol;

Now we begin the nonlinear loop

                for (unsigned int l=0; l<n_nonlinear_steps; ++l)
        	  {

Update the nonlinear solution.

                    last_nonlinear_soln->zero();
        	    last_nonlinear_soln->add(*navier_stokes_system.solution);

Assemble & solve the linear system.

                    perf_log.start_event("linear solve");
        	    equation_systems.get_system("Navier-Stokes").solve();
        	    perf_log.stop_event("linear solve");

Compute the difference between this solution and the last nonlinear iterate.

                    last_nonlinear_soln->add (-1., *navier_stokes_system.solution);

Close the vector before computing its norm

                    last_nonlinear_soln->close();

Compute the l2 norm of the difference

                    const Real norm_delta = last_nonlinear_soln->l2_norm();

How many iterations were required to solve the linear system?

                    const unsigned int n_linear_iterations = navier_stokes_system.n_linear_iterations();

What was the final residual of the linear system?

                    const Real final_linear_residual = navier_stokes_system.final_linear_residual();

Print out convergence information for the linear and nonlinear iterations.

                    std::cout << "Linear solver converged at step: "
        		      << n_linear_iterations
        		      << ", final residual: "
        		      << final_linear_residual
        		      << "  Nonlinear convergence: ||u - u_old|| = "
        		      << norm_delta
        		      << std::endl;

Terminate the solution iteration if the difference between this nonlinear iterate and the last is sufficiently small, AND if the most recent linear system was solved to a sufficient tolerance.

                    if ((norm_delta < nonlinear_tolerance) &&
        		(navier_stokes_system.final_linear_residual() < nonlinear_tolerance))
        	      {
        		std::cout << " Nonlinear solver converged at step "
        			  << l
        			  << std::endl;
        		break;
        	      }

Otherwise, decrease the linear system tolerance. For the inexact Newton method, the linear solver tolerance needs to decrease as we get closer to the solution to ensure quadratic convergence. The new linear solver tolerance is chosen (heuristically) as the square of the previous linear system residual norm. Real flr2 = final_linear_residual*final_linear_residual;

                    equation_systems.parameters.set<Real> ("linear solver tolerance") =
        	      Utility::pow<2>(final_linear_residual);
        
        	  } // end nonlinear loop

Write out every nth timestep to file.

                const unsigned int write_interval = 1;
        	
        	if ((t_step+1)%write_interval == 0)
        	  {
        	    OStringStream file_name;

We write the file name in the gmv auto-read format.

                    file_name << "out.gmv.";
        	    OSSRealzeroright(file_name,3,0, t_step + 1);
        	    
        	    GMVIO(mesh).write_equation_systems (file_name.str(),
        						equation_systems);
        	  }
              } // end timestep loop.
          }

All done.

          return libMesh::close ();
        }

The matrix assembly function to be called at each time step to prepare for the linear solve.

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

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

          assert (system_name == "Navier-Stokes");

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 Stokes system object.

          TransientLinearImplicitSystem & navier_stokes_system =
            es.get_system<TransientLinearImplicitSystem> ("Navier-Stokes");

Numeric ids corresponding to each variable in the system

          const unsigned int u_var = navier_stokes_system.variable_number ("u");
          const unsigned int v_var = navier_stokes_system.variable_number ("v");
          const unsigned int p_var = navier_stokes_system.variable_number ("p");

Get the Finite Element type for "u". Note this will be the same as the type for "v".

          FEType fe_vel_type = navier_stokes_system.variable_type(u_var);

Get the Finite Element type for "p".

          FEType fe_pres_type = navier_stokes_system.variable_type(p_var);

Build a Finite Element object of the specified type for the velocity variables.

          AutoPtr<FEBase> fe_vel  (FEBase::build(dim, fe_vel_type));

Build a Finite Element object of the specified type for the pressure variables.

          AutoPtr<FEBase> fe_pres (FEBase::build(dim, fe_pres_type));

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

          QGauss qrule (dim, fe_vel_type.default_quadrature_order());

Tell the finite element objects to use our quadrature rule.

          fe_vel->attach_quadrature_rule (&qrule);
          fe_pres->attach_quadrature_rule (&qrule);

Here we define some references to cell-specific data that will be used to assemble the linear system.

The element Jacobian * quadrature weight at each integration point.

          const std::vector<Real>& JxW = fe_vel->get_JxW();

The element shape functions evaluated at the quadrature points.

          const std::vector<std::vector<Real> >& phi = fe_vel->get_phi();

The element shape function gradients for the velocity variables evaluated at the quadrature points.

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

The element shape functions for the pressure variable evaluated at the quadrature points.

          const std::vector<std::vector<Real> >& psi = fe_pres->get_phi();

The value of the linear shape function gradients at the quadrature points const std::vector >& dpsi = fe_pres->get_dphi();

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 = navier_stokes_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;
        
          DenseSubMatrix<Number>
            Kuu(Ke), Kuv(Ke), Kup(Ke),
            Kvu(Ke), Kvv(Ke), Kvp(Ke),
            Kpu(Ke), Kpv(Ke), Kpp(Ke);
        
          DenseSubVector<Number>
            Fu(Fe),
            Fv(Fe),
            Fp(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;
          std::vector<unsigned int> dof_indices_u;
          std::vector<unsigned int> dof_indices_v;
          std::vector<unsigned int> dof_indices_p;

Find out what the timestep size parameter is from the system, and the value of theta for the theta method. We use implicit Euler (theta=1) for this simulation even though it is only first-order accurate in time. The reason for this decision is that the second-order Crank-Nicolson method is notoriously oscillatory for problems with discontinuous initial data such as the lid-driven cavity. Therefore, we sacrifice accuracy in time for stability, but since the solution reaches steady state relatively quickly we can afford to take small timesteps. If you monitor the initial nonlinear residual for this simulation, you should see that it is monotonically decreasing in time.

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

const Real time = es.parameters.get("time");

          const Real theta = 1.;

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.

          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);
              dof_map.dof_indices (elem, dof_indices_u, u_var);
              dof_map.dof_indices (elem, dof_indices_v, v_var);
              dof_map.dof_indices (elem, dof_indices_p, p_var);
        
              const unsigned int n_dofs   = dof_indices.size();
              const unsigned int n_u_dofs = dof_indices_u.size(); 
              const unsigned int n_v_dofs = dof_indices_v.size(); 
              const unsigned int n_p_dofs = dof_indices_p.size();

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_vel->reinit  (elem);
              fe_pres->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 (n_dofs, n_dofs);
              Fe.resize (n_dofs);

Reposition the submatrices... The idea is this:

- - - - | Kuu Kuv Kup | | Fu | Ke = | Kvu Kvv Kvp |; Fe = | Fv | | Kpu Kpv Kpp | | Fp | - - - -

The \p DenseSubMatrix.repostition () member takes the (row_offset, column_offset, row_size, column_size).

Similarly, the \p DenseSubVector.reposition () member takes the (row_offset, row_size)

              Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs);
              Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs);
              Kup.reposition (u_var*n_u_dofs, p_var*n_u_dofs, n_u_dofs, n_p_dofs);
              
              Kvu.reposition (v_var*n_v_dofs, u_var*n_v_dofs, n_v_dofs, n_u_dofs);
              Kvv.reposition (v_var*n_v_dofs, v_var*n_v_dofs, n_v_dofs, n_v_dofs);
              Kvp.reposition (v_var*n_v_dofs, p_var*n_v_dofs, n_v_dofs, n_p_dofs);
              
              Kpu.reposition (p_var*n_u_dofs, u_var*n_u_dofs, n_p_dofs, n_u_dofs);
              Kpv.reposition (p_var*n_u_dofs, v_var*n_u_dofs, n_p_dofs, n_v_dofs);
              Kpp.reposition (p_var*n_u_dofs, p_var*n_u_dofs, n_p_dofs, n_p_dofs);
        
              Fu.reposition (u_var*n_u_dofs, n_u_dofs);
              Fv.reposition (v_var*n_u_dofs, n_v_dofs);
              Fp.reposition (p_var*n_u_dofs, n_p_dofs);

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 must 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 solution & its gradient at the previous timestep.

                  Number   u = 0., u_old = 0.;
        	  Number   v = 0., v_old = 0.;
        	  Number   p_old = 0.;
        	  Gradient grad_u, grad_u_old;
        	  Gradient grad_v, grad_v_old;

Compute the velocity & its gradient from the previous timestep and the old Newton iterate.

                  for (unsigned int l=0; l<n_u_dofs; l++)
        	    {

From the old timestep:

                      u_old += phi[l][qp]*navier_stokes_system.old_solution (dof_indices_u[l]);
        	      v_old += phi[l][qp]*navier_stokes_system.old_solution (dof_indices_v[l]);
        	      grad_u_old.add_scaled (dphi[l][qp],navier_stokes_system.old_solution (dof_indices_u[l]));
        	      grad_v_old.add_scaled (dphi[l][qp],navier_stokes_system.old_solution (dof_indices_v[l]));

From the previous Newton iterate:

                      u += phi[l][qp]*navier_stokes_system.current_solution (dof_indices_u[l]); 
        	      v += phi[l][qp]*navier_stokes_system.current_solution (dof_indices_v[l]);
        	      grad_u.add_scaled (dphi[l][qp],navier_stokes_system.current_solution (dof_indices_u[l]));
        	      grad_v.add_scaled (dphi[l][qp],navier_stokes_system.current_solution (dof_indices_v[l]));
        	    }

Compute the old pressure value at this quadrature point.

                  for (unsigned int l=0; l<n_p_dofs; l++)
        	    {
        	      p_old += psi[l][qp]*navier_stokes_system.old_solution (dof_indices_p[l]);
        	    }

Definitions for convenience. It is sometimes simpler to do a dot product if you have the full vector at your disposal.

                  const NumberVectorValue U_old (u_old, v_old);
        	  const NumberVectorValue U     (u,     v);
        	  const Number  u_x = grad_u(0);
        	  const Number  u_y = grad_u(1);
        	  const Number  v_x = grad_v(0);
        	  const Number  v_y = grad_v(1);

First, an i-loop over the velocity degrees of freedom. We know that n_u_dofs == n_v_dofs so we can compute contributions for both at the same time.

                  for (unsigned int i=0; i<n_u_dofs; i++)
        	    {
        	      Fu(i) += JxW[qp]*(u_old*phi[i][qp] -                            // mass-matrix term 
        				(1.-theta)*dt*(U_old*grad_u_old)*phi[i][qp] + // convection term
        				(1.-theta)*dt*p_old*dphi[i][qp](0)  -         // pressure term on rhs
        				(1.-theta)*dt*(grad_u_old*dphi[i][qp]) +      // diffusion term on rhs
        				theta*dt*(U*grad_u)*phi[i][qp]);              // Newton term
        
        		
        	      Fv(i) += JxW[qp]*(v_old*phi[i][qp] -                             // mass-matrix term
        				(1.-theta)*dt*(U_old*grad_v_old)*phi[i][qp] +  // convection term
        				(1.-theta)*dt*p_old*dphi[i][qp](1) -           // pressure term on rhs
        				(1.-theta)*dt*(grad_v_old*dphi[i][qp]) +       // diffusion term on rhs
        				theta*dt*(U*grad_v)*phi[i][qp]);               // Newton term

Note that the Fp block is identically zero unless we are using some kind of artificial compressibility scheme...

Matrix contributions for the uu and vv couplings.

                      for (unsigned int j=0; j<n_u_dofs; j++)
        		{
        		  Kuu(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp] +                // mass matrix term
        				       theta*dt*(dphi[i][qp]*dphi[j][qp]) +   // diffusion term
        				       theta*dt*(U*dphi[j][qp])*phi[i][qp] +  // convection term
        				       theta*dt*u_x*phi[i][qp]*phi[j][qp]);   // Newton term
        
        		  Kuv(i,j) += JxW[qp]*theta*dt*u_y*phi[i][qp]*phi[j][qp];     // Newton term
        		  
        		  Kvv(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp] +                // mass matrix term
        				       theta*dt*(dphi[i][qp]*dphi[j][qp]) +   // diffusion term
        				       theta*dt*(U*dphi[j][qp])*phi[i][qp] +  // convection term
        				       theta*dt*v_y*phi[i][qp]*phi[j][qp]);   // Newton term
        
        		  Kvu(i,j) += JxW[qp]*theta*dt*v_x*phi[i][qp]*phi[j][qp];     // Newton term
        		}

Matrix contributions for the up and vp couplings.

                      for (unsigned int j=0; j<n_p_dofs; j++)
        		{
        		  Kup(i,j) += JxW[qp]*(-theta*dt*psi[j][qp]*dphi[i][qp](0));
        		  Kvp(i,j) += JxW[qp]*(-theta*dt*psi[j][qp]*dphi[i][qp](1));
        		}
        	    }

Now an i-loop over the pressure degrees of freedom. This code computes the matrix entries due to the continuity equation. Note: To maintain a symmetric matrix, we may (or may not) multiply the continuity equation by negative one. Here we do not.

                  for (unsigned int i=0; i<n_p_dofs; i++)
        	    {
        	      for (unsigned int j=0; j<n_u_dofs; j++)
        		{
        		  Kpu(i,j) += JxW[qp]*psi[i][qp]*dphi[j][qp](0);
        		  Kpv(i,j) += JxW[qp]*psi[i][qp]*dphi[j][qp](1);
        		}
        	    }
        	} // end of the quadrature point qp-loop

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 penalty method used here is equivalent (for Lagrange basis functions) to lumping the matrix resulting from the L2 projection penalty approach introduced in example 3.

The penalty value. \f$ \frac{1}{\epsilon \f$

                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)
        	    {

Get the boundary ID for side 's'. These are set internally by build_square(). 0=bottom 1=right 2=top 3=left

                      short int bc_id = mesh.boundary_info->boundary_id (elem,s);
        	      if (bc_id==BoundaryInfo::invalid_id)
        		  error();
        
        	      
        	      AutoPtr<Elem> side (elem->build_side(s));

Loop over the nodes on the side.

                      for (unsigned int ns=0; ns<side->n_nodes(); ns++)
        		{

Get the boundary values.

Set u = 1 on the top boundary, 0 everywhere else

                          const Real u_value = (bc_id==2) ? 1. : 0.;

Set v = 0 everywhere

                          const Real v_value = 0.;

Find the node on the element matching this node on the side. That defined where in the element matrix the boundary condition will be applied.

                          for (unsigned int n=0; n<elem->n_nodes(); n++)
        		    if (elem->node(n) == side->node(ns))
        		      {

Matrix contribution.

                                Kuu(n,n) += penalty;
        			Kvv(n,n) += penalty;

Right-hand-side contribution.

                                Fu(n) += penalty*u_value;
        			Fv(n) += penalty*v_value;
        		      }
        		} // end face node loop	  
        	    } // end if (elem->neighbor(side) == NULL)

Pin the pressure to zero at global node number "pressure_node". This effectively removes the non-trivial null space of constant pressure solutions.

                const bool pin_pressure = true;
        	if (pin_pressure)
        	  {
        	    const unsigned int pressure_node = 0;
        	    const Real p_value               = 0.0;
        	    for (unsigned int c=0; c<elem->n_nodes(); c++)
        	      if (elem->node(c) == pressure_node)
        		{
        		  Kpp(c,c) += penalty;
        		  Fp(c)    += penalty*p_value;
        		}
        	  }
              } // end boundary condition section

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.

              navier_stokes_system.matrix->add_matrix (Ke, dof_indices);
              navier_stokes_system.rhs->add_vector    (Fe, dof_indices);
            } // end of element loop

That's it.

          return;
        }

The program without comments:

 
  
  #include <iostream>
  #include <algorithm>
  #include <math.h>
  
  #include "libmesh.h"
  #include "mesh.h"
  #include "mesh_generation.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 "linear_implicit_system.h"
  #include "transient_system.h"
  #include "perf_log.h"
  #include "boundary_info.h"
  #include "utility.h"
  
  #include "o_string_stream.h"
  
  #include "dense_submatrix.h"
  #include "dense_subvector.h"
  
  #include "elem.h"
  
  void assemble_stokes (EquationSystems& es,
  		      const std::string& system_name);
  
  int main (int argc, char** argv)
  {
    libMesh::init (argc, argv);
    {    
      const unsigned int dim = 2;     
      
      Mesh mesh (dim);
      
      MeshTools::Generation::build_square (mesh,
  					 20, 20,
  					 0., 1.,
  					 0., 1.,
  					 QUAD9);
      
      mesh.print_info();
      
      EquationSystems equation_systems (mesh);
      
      {
        TransientLinearImplicitSystem & system = 
  	equation_systems.add_system<TransientLinearImplicitSystem> ("Navier-Stokes");
        
        system.add_variable ("u", SECOND);
        system.add_variable ("v", SECOND);
  
        system.add_variable ("p", FIRST);
  
        system.attach_assemble_function (assemble_stokes);
        
        equation_systems.init ();
  
        equation_systems.print_info();
      }
  
      PerfLog perf_log("Example 13");
      
      const Real dt = 0.01;
      Real time     = 0.0;
      const unsigned int n_timesteps = 15;
  
      const unsigned int n_nonlinear_steps = 15;
      const Real nonlinear_tolerance       = 1.e-3;
  
      equation_systems.parameters.set<unsigned int>("linear solver maximum iterations") = 250;
      
      equation_systems.parameters.set<Real> ("dt")   = dt;
  
      TransientLinearImplicitSystem&  navier_stokes_system =
  	  equation_systems.get_system<TransientLinearImplicitSystem>("Navier-Stokes");
  
      AutoPtr<NumericVector<Number> >
        last_nonlinear_soln (navier_stokes_system.solution->clone());
  
      for (unsigned int t_step=0; t_step<n_timesteps; ++t_step)
        {
  	time += dt;
  
  	equation_systems.parameters.set<Real> ("time") = time;
  
  	std::cout << "\n\n*** Solving time step " << t_step << ", time = " << time << " ***" << std::endl;
  
  	*navier_stokes_system.old_local_solution = *navier_stokes_system.current_local_solution;
  
  	const Real initial_linear_solver_tol = 1.e-6;
  	equation_systems.parameters.set<Real> ("linear solver tolerance") = initial_linear_solver_tol;
  
  	for (unsigned int l=0; l<n_nonlinear_steps; ++l)
  	  {
  	    last_nonlinear_soln->zero();
  	    last_nonlinear_soln->add(*navier_stokes_system.solution);
  	    
  	    perf_log.start_event("linear solve");
  	    equation_systems.get_system("Navier-Stokes").solve();
  	    perf_log.stop_event("linear solve");
  
  	    last_nonlinear_soln->add (-1., *navier_stokes_system.solution);
  
  	    last_nonlinear_soln->close();
  
  	    const Real norm_delta = last_nonlinear_soln->l2_norm();
  
  	    const unsigned int n_linear_iterations = navier_stokes_system.n_linear_iterations();
  	    
  	    const Real final_linear_residual = navier_stokes_system.final_linear_residual();
  	    
  	    std::cout << "Linear solver converged at step: "
  		      << n_linear_iterations
  		      << ", final residual: "
  		      << final_linear_residual
  		      << "  Nonlinear convergence: ||u - u_old|| = "
  		      << norm_delta
  		      << std::endl;
  
  	    if ((norm_delta < nonlinear_tolerance) &&
  		(navier_stokes_system.final_linear_residual() < nonlinear_tolerance))
  	      {
  		std::cout << " Nonlinear solver converged at step "
  			  << l
  			  << std::endl;
  		break;
  	      }
  	    
  	    equation_systems.parameters.set<Real> ("linear solver tolerance") =
  	      Utility::pow<2>(final_linear_residual);
  
  	  } // end nonlinear loop
  	
  	const unsigned int write_interval = 1;
  	
  	if ((t_step+1)%write_interval == 0)
  	  {
  	    OStringStream file_name;
  
  	    file_name << "out.gmv.";
  	    OSSRealzeroright(file_name,3,0, t_step + 1);
  	    
  	    GMVIO(mesh).write_equation_systems (file_name.str(),
  						equation_systems);
  	  }
        } // end timestep loop.
    }
    
    return libMesh::close ();
  }
  
  
  
  
  
  
  void assemble_stokes (EquationSystems& es,
  		      const std::string& system_name)
  {
    assert (system_name == "Navier-Stokes");
    
    const Mesh& mesh = es.get_mesh();
    
    const unsigned int dim = mesh.mesh_dimension();
    
    TransientLinearImplicitSystem & navier_stokes_system =
      es.get_system<TransientLinearImplicitSystem> ("Navier-Stokes");
  
    const unsigned int u_var = navier_stokes_system.variable_number ("u");
    const unsigned int v_var = navier_stokes_system.variable_number ("v");
    const unsigned int p_var = navier_stokes_system.variable_number ("p");
    
    FEType fe_vel_type = navier_stokes_system.variable_type(u_var);
    
    FEType fe_pres_type = navier_stokes_system.variable_type(p_var);
  
    AutoPtr<FEBase> fe_vel  (FEBase::build(dim, fe_vel_type));
      
    AutoPtr<FEBase> fe_pres (FEBase::build(dim, fe_pres_type));
    
    QGauss qrule (dim, fe_vel_type.default_quadrature_order());
  
    fe_vel->attach_quadrature_rule (&qrule);
    fe_pres->attach_quadrature_rule (&qrule);
    
    const std::vector<Real>& JxW = fe_vel->get_JxW();
  
    const std::vector<std::vector<Real> >& phi = fe_vel->get_phi();
  
    const std::vector<std::vector<RealGradient> >& dphi = fe_vel->get_dphi();
  
    const std::vector<std::vector<Real> >& psi = fe_pres->get_phi();
  
    
    const DofMap & dof_map = navier_stokes_system.get_dof_map();
  
    DenseMatrix<Number> Ke;
    DenseVector<Number> Fe;
  
    DenseSubMatrix<Number>
      Kuu(Ke), Kuv(Ke), Kup(Ke),
      Kvu(Ke), Kvv(Ke), Kvp(Ke),
      Kpu(Ke), Kpv(Ke), Kpp(Ke);
  
    DenseSubVector<Number>
      Fu(Fe),
      Fv(Fe),
      Fp(Fe);
  
    std::vector<unsigned int> dof_indices;
    std::vector<unsigned int> dof_indices_u;
    std::vector<unsigned int> dof_indices_v;
    std::vector<unsigned int> dof_indices_p;
  
    const Real dt    = es.parameters.get<Real>("dt");
    const Real theta = 1.;
      
    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);
        dof_map.dof_indices (elem, dof_indices_u, u_var);
        dof_map.dof_indices (elem, dof_indices_v, v_var);
        dof_map.dof_indices (elem, dof_indices_p, p_var);
  
        const unsigned int n_dofs   = dof_indices.size();
        const unsigned int n_u_dofs = dof_indices_u.size(); 
        const unsigned int n_v_dofs = dof_indices_v.size(); 
        const unsigned int n_p_dofs = dof_indices_p.size();
        
        fe_vel->reinit  (elem);
        fe_pres->reinit (elem);
  
        Ke.resize (n_dofs, n_dofs);
        Fe.resize (n_dofs);
  
        Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs);
        Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs);
        Kup.reposition (u_var*n_u_dofs, p_var*n_u_dofs, n_u_dofs, n_p_dofs);
        
        Kvu.reposition (v_var*n_v_dofs, u_var*n_v_dofs, n_v_dofs, n_u_dofs);
        Kvv.reposition (v_var*n_v_dofs, v_var*n_v_dofs, n_v_dofs, n_v_dofs);
        Kvp.reposition (v_var*n_v_dofs, p_var*n_v_dofs, n_v_dofs, n_p_dofs);
        
        Kpu.reposition (p_var*n_u_dofs, u_var*n_u_dofs, n_p_dofs, n_u_dofs);
        Kpv.reposition (p_var*n_u_dofs, v_var*n_u_dofs, n_p_dofs, n_v_dofs);
        Kpp.reposition (p_var*n_u_dofs, p_var*n_u_dofs, n_p_dofs, n_p_dofs);
  
        Fu.reposition (u_var*n_u_dofs, n_u_dofs);
        Fv.reposition (v_var*n_u_dofs, n_v_dofs);
        Fp.reposition (p_var*n_u_dofs, n_p_dofs);
  
        for (unsigned int qp=0; qp<qrule.n_points(); qp++)
  	{
  	  Number   u = 0., u_old = 0.;
  	  Number   v = 0., v_old = 0.;
  	  Number   p_old = 0.;
  	  Gradient grad_u, grad_u_old;
  	  Gradient grad_v, grad_v_old;
  	  
  	  for (unsigned int l=0; l<n_u_dofs; l++)
  	    {
  	      u_old += phi[l][qp]*navier_stokes_system.old_solution (dof_indices_u[l]);
  	      v_old += phi[l][qp]*navier_stokes_system.old_solution (dof_indices_v[l]);
  	      grad_u_old.add_scaled (dphi[l][qp],navier_stokes_system.old_solution (dof_indices_u[l]));
  	      grad_v_old.add_scaled (dphi[l][qp],navier_stokes_system.old_solution (dof_indices_v[l]));
  
  	      u += phi[l][qp]*navier_stokes_system.current_solution (dof_indices_u[l]); 
  	      v += phi[l][qp]*navier_stokes_system.current_solution (dof_indices_v[l]);
  	      grad_u.add_scaled (dphi[l][qp],navier_stokes_system.current_solution (dof_indices_u[l]));
  	      grad_v.add_scaled (dphi[l][qp],navier_stokes_system.current_solution (dof_indices_v[l]));
  	    }
  
  	  for (unsigned int l=0; l<n_p_dofs; l++)
  	    {
  	      p_old += psi[l][qp]*navier_stokes_system.old_solution (dof_indices_p[l]);
  	    }
  
  	  const NumberVectorValue U_old (u_old, v_old);
  	  const NumberVectorValue U     (u,     v);
  	  const Number  u_x = grad_u(0);
  	  const Number  u_y = grad_u(1);
  	  const Number  v_x = grad_v(0);
  	  const Number  v_y = grad_v(1);
  	  
  	  for (unsigned int i=0; i<n_u_dofs; i++)
  	    {
  	      Fu(i) += JxW[qp]*(u_old*phi[i][qp] -                            // mass-matrix term 
  				(1.-theta)*dt*(U_old*grad_u_old)*phi[i][qp] + // convection term
  				(1.-theta)*dt*p_old*dphi[i][qp](0)  -         // pressure term on rhs
  				(1.-theta)*dt*(grad_u_old*dphi[i][qp]) +      // diffusion term on rhs
  				theta*dt*(U*grad_u)*phi[i][qp]);              // Newton term
  
  		
  	      Fv(i) += JxW[qp]*(v_old*phi[i][qp] -                             // mass-matrix term
  				(1.-theta)*dt*(U_old*grad_v_old)*phi[i][qp] +  // convection term
  				(1.-theta)*dt*p_old*dphi[i][qp](1) -           // pressure term on rhs
  				(1.-theta)*dt*(grad_v_old*dphi[i][qp]) +       // diffusion term on rhs
  				theta*dt*(U*grad_v)*phi[i][qp]);               // Newton term
  	    
  
  
  	      for (unsigned int j=0; j<n_u_dofs; j++)
  		{
  		  Kuu(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp] +                // mass matrix term
  				       theta*dt*(dphi[i][qp]*dphi[j][qp]) +   // diffusion term
  				       theta*dt*(U*dphi[j][qp])*phi[i][qp] +  // convection term
  				       theta*dt*u_x*phi[i][qp]*phi[j][qp]);   // Newton term
  
  		  Kuv(i,j) += JxW[qp]*theta*dt*u_y*phi[i][qp]*phi[j][qp];     // Newton term
  		  
  		  Kvv(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp] +                // mass matrix term
  				       theta*dt*(dphi[i][qp]*dphi[j][qp]) +   // diffusion term
  				       theta*dt*(U*dphi[j][qp])*phi[i][qp] +  // convection term
  				       theta*dt*v_y*phi[i][qp]*phi[j][qp]);   // Newton term
  
  		  Kvu(i,j) += JxW[qp]*theta*dt*v_x*phi[i][qp]*phi[j][qp];     // Newton term
  		}
  
  	      for (unsigned int j=0; j<n_p_dofs; j++)
  		{
  		  Kup(i,j) += JxW[qp]*(-theta*dt*psi[j][qp]*dphi[i][qp](0));
  		  Kvp(i,j) += JxW[qp]*(-theta*dt*psi[j][qp]*dphi[i][qp](1));
  		}
  	    }
  
  	  for (unsigned int i=0; i<n_p_dofs; i++)
  	    {
  	      for (unsigned int j=0; j<n_u_dofs; j++)
  		{
  		  Kpu(i,j) += JxW[qp]*psi[i][qp]*dphi[j][qp](0);
  		  Kpv(i,j) += JxW[qp]*psi[i][qp]*dphi[j][qp](1);
  		}
  	    }
  	} // end of the quadrature point qp-loop
  
        
        {
  	const Real penalty = 1.e10;
  		  
  	for (unsigned int s=0; s<elem->n_sides(); s++)
  	  if (elem->neighbor(s) == NULL)
  	    {
  	      short int bc_id = mesh.boundary_info->boundary_id (elem,s);
  	      if (bc_id==BoundaryInfo::invalid_id)
  		  error();
  
  	      
  	      AutoPtr<Elem> side (elem->build_side(s));
  	      	      
  	      for (unsigned int ns=0; ns<side->n_nodes(); ns++)
  		{
  		   
  		  const Real u_value = (bc_id==2) ? 1. : 0.;
  		  
  		  const Real v_value = 0.;
  		  
  		  for (unsigned int n=0; n<elem->n_nodes(); n++)
  		    if (elem->node(n) == side->node(ns))
  		      {
  			Kuu(n,n) += penalty;
  			Kvv(n,n) += penalty;
  		  		  
  			Fu(n) += penalty*u_value;
  			Fv(n) += penalty*v_value;
  		      }
  		} // end face node loop	  
  	    } // end if (elem->neighbor(side) == NULL)
  	
  	const bool pin_pressure = true;
  	if (pin_pressure)
  	  {
  	    const unsigned int pressure_node = 0;
  	    const Real p_value               = 0.0;
  	    for (unsigned int c=0; c<elem->n_nodes(); c++)
  	      if (elem->node(c) == pressure_node)
  		{
  		  Kpp(c,c) += penalty;
  		  Fp(c)    += penalty*p_value;
  		}
  	  }
        } // end boundary condition section	  
        
        navier_stokes_system.matrix->add_matrix (Ke, dof_indices);
        navier_stokes_system.rhs->add_vector    (Fe, dof_indices);
      } // end of element loop
    
    return;
  }

The console output of the program:

***************************************************************
* Running Example  ./ex13-devel
***************************************************************
 
 Mesh Information:
  mesh_dimension()=2
  spatial_dimension()=3
  n_nodes()=1681
  n_elem()=400
   n_local_elem()=400
   n_active_elem()=400
  n_subdomains()=1
  n_processors()=1
  processor_id()=0

 EquationSystems
  n_systems()=1
   System "Navier-Stokes"
    Type "TransientLinearImplicit"
    Variables="u" "v" "p" 
    Finite Element Types="LAGRANGE" "LAGRANGE" "LAGRANGE" 
    Approximation Orders="SECOND" "SECOND" "FIRST" 
    n_dofs()=3803
    n_local_dofs()=3803
    n_constrained_dofs()=0
    n_vectors()=3



*** Solving time step 0, time = 0.01 ***
Linear solver converged at step: 236, final residual: 0.000214866  Nonlinear convergence: ||u - u_old|| = 282.453
Linear solver converged at step: 90, final residual: 9.99752e-06  Nonlinear convergence: ||u - u_old|| = 0.743588
Linear solver converged at step: 250, final residual: 1.04334e-07  Nonlinear convergence: ||u - u_old|| = 0.0151695
Linear solver converged at step: 250, final residual: 4.20802e-10  Nonlinear convergence: ||u - u_old|| = 0.000178465
 Nonlinear solver converged at step 3


*** Solving time step 1, time = 0.02 ***
Linear solver converged at step: 148, final residual: 0.000196443  Nonlinear convergence: ||u - u_old|| = 30.6968
Linear solver converged at step: 74, final residual: 8.13638e-06  Nonlinear convergence: ||u - u_old|| = 0.194999
Linear solver converged at step: 250, final residual: 3.57667e-08  Nonlinear convergence: ||u - u_old|| = 0.0088697
Linear solver converged at step: 250, final residual: 8.44767e-13  Nonlinear convergence: ||u - u_old|| = 5.9269e-05
 Nonlinear solver converged at step 3


*** Solving time step 2, time = 0.03 ***
Linear solver converged at step: 77, final residual: 0.000223335  Nonlinear convergence: ||u - u_old|| = 6.00665
Linear solver converged at step: 118, final residual: 9.30581e-06  Nonlinear convergence: ||u - u_old|| = 0.268524
Linear solver converged at step: 250, final residual: 3.29197e-08  Nonlinear convergence: ||u - u_old|| = 0.00316201
Linear solver converged at step: 250, final residual: 2.45713e-09  Nonlinear convergence: ||u - u_old|| = 5.82116e-05
 Nonlinear solver converged at step 3


*** Solving time step 3, time = 0.04 ***
Linear solver converged at step: 72, final residual: 0.000210495  Nonlinear convergence: ||u - u_old|| = 2.18249
Linear solver converged at step: 159, final residual: 9.538e-06  Nonlinear convergence: ||u - u_old|| = 0.260483
Linear solver converged at step: 250, final residual: 2.18784e-07  Nonlinear convergence: ||u - u_old|| = 0.0103119
Linear solver converged at step: 250, final residual: 3.4177e-09  Nonlinear convergence: ||u - u_old|| = 0.000386804
 Nonlinear solver converged at step 3


*** Solving time step 4, time = 0.05 ***
Linear solver converged at step: 58, final residual: 0.000192962  Nonlinear convergence: ||u - u_old|| = 1.02052
Linear solver converged at step: 119, final residual: 7.33031e-06  Nonlinear convergence: ||u - u_old|| = 0.103956
Linear solver converged at step: 225, final residual: 1.16591e-08  Nonlinear convergence: ||u - u_old|| = 0.00534628
Linear solver converged at step: 250, final residual: 1.41694e-12  Nonlinear convergence: ||u - u_old|| = 3.60742e-06
 Nonlinear solver converged at step 3


*** Solving time step 5, time = 0.06 ***
Linear solver converged at step: 52, final residual: 0.000214167  Nonlinear convergence: ||u - u_old|| = 0.476324
Linear solver converged at step: 172, final residual: 1.00119e-05  Nonlinear convergence: ||u - u_old|| = 0.202403
Linear solver converged at step: 158, final residual: 2.16565e-08  Nonlinear convergence: ||u - u_old|| = 0.0173926
Linear solver converged at step: 250, final residual: 2.6449e-11  Nonlinear convergence: ||u - u_old|| = 2.11892e-05
 Nonlinear solver converged at step 3


*** Solving time step 6, time = 0.07 ***
Linear solver converged at step: 21, final residual: 0.000225026  Nonlinear convergence: ||u - u_old|| = 0.256854
Linear solver converged at step: 117, final residual: 1.12946e-05  Nonlinear convergence: ||u - u_old|| = 0.147309
Linear solver converged at step: 250, final residual: 1.54091e-06  Nonlinear convergence: ||u - u_old|| = 0.0128016
Linear solver converged at step: 250, final residual: 2.11486e-09  Nonlinear convergence: ||u - u_old|| = 0.00239393
Linear solver converged at step: 250, final residual: 8.06203e-12  Nonlinear convergence: ||u - u_old|| = 2.67758e-06
 Nonlinear solver converged at step 4


*** Solving time step 7, time = 0.08 ***
Linear solver converged at step: 19, final residual: 0.000192228  Nonlinear convergence: ||u - u_old|| = 0.152223
Linear solver converged at step: 67, final residual: 7.02725e-06  Nonlinear convergence: ||u - u_old|| = 0.0850415
Linear solver converged at step: 194, final residual: 1.08592e-08  Nonlinear convergence: ||u - u_old|| = 0.00357155
Linear solver converged at step: 250, final residual: 6.3798e-12  Nonlinear convergence: ||u - u_old|| = 1.67316e-05
 Nonlinear solver converged at step 3


*** Solving time step 8, time = 0.09 ***
Linear solver converged at step: 15, final residual: 0.000204876  Nonlinear convergence: ||u - u_old|| = 0.086678
Linear solver converged at step: 97, final residual: 8.76745e-06  Nonlinear convergence: ||u - u_old|| = 0.0454467
Linear solver converged at step: 137, final residual: 1.72036e-08  Nonlinear convergence: ||u - u_old|| = 0.00314358
Linear solver converged at step: 250, final residual: 3.37166e-12  Nonlinear convergence: ||u - u_old|| = 1.96889e-05
 Nonlinear solver converged at step 3


*** Solving time step 9, time = 0.1 ***
Linear solver converged at step: 11, final residual: 0.00022001  Nonlinear convergence: ||u - u_old|| = 0.0505868
Linear solver converged at step: 69, final residual: 1.08054e-05  Nonlinear convergence: ||u - u_old|| = 0.0182682
Linear solver converged at step: 250, final residual: 3.49658e-08  Nonlinear convergence: ||u - u_old|| = 0.0127024
Linear solver converged at step: 250, final residual: 2.95441e-10  Nonlinear convergence: ||u - u_old|| = 5.53751e-05
 Nonlinear solver converged at step 3


*** Solving time step 10, time = 0.11 ***
Linear solver converged at step: 11, final residual: 0.000139574  Nonlinear convergence: ||u - u_old|| = 0.0322373
Linear solver converged at step: 104, final residual: 4.33263e-06  Nonlinear convergence: ||u - u_old|| = 0.0126853
Linear solver converged at step: 250, final residual: 4.09989e-08  Nonlinear convergence: ||u - u_old|| = 0.00565014
Linear solver converged at step: 250, final residual: 1.93796e-10  Nonlinear convergence: ||u - u_old|| = 5.64589e-05
 Nonlinear solver converged at step 3


*** Solving time step 11, time = 0.12 ***
Linear solver converged at step: 10, final residual: 0.000145884  Nonlinear convergence: ||u - u_old|| = 0.0205125
Linear solver converged at step: 48, final residual: 4.75062e-06  Nonlinear convergence: ||u - u_old|| = 0.00781029
Linear solver converged at step: 250, final residual: 1.24408e-08  Nonlinear convergence: ||u - u_old|| = 0.00689828
Linear solver converged at step: 250, final residual: 1.17963e-10  Nonlinear convergence: ||u - u_old|| = 2.01527e-05
 Nonlinear solver converged at step 3


*** Solving time step 12, time = 0.13 ***
Linear solver converged at step: 10, final residual: 9.34716e-05  Nonlinear convergence: ||u - u_old|| = 0.0136016
Linear solver converged at step: 68, final residual: 1.95668e-06  Nonlinear convergence: ||u - u_old|| = 0.00578812
Linear solver converged at step: 250, final residual: 4.84151e-09  Nonlinear convergence: ||u - u_old|| = 0.00253899
Linear solver converged at step: 250, final residual: 3.95845e-11  Nonlinear convergence: ||u - u_old|| = 9.45948e-06
 Nonlinear solver converged at step 3


*** Solving time step 13, time = 0.14 ***
Linear solver converged at step: 9, final residual: 0.000149994  Nonlinear convergence: ||u - u_old|| = 0.00910172
Linear solver converged at step: 19, final residual: 3.88794e-06  Nonlinear convergence: ||u - u_old|| = 0.00478718
Linear solver converged at step: 250, final residual: 1.09364e-08  Nonlinear convergence: ||u - u_old|| = 0.00216678
Linear solver converged at step: 250, final residual: 5.72984e-10  Nonlinear convergence: ||u - u_old|| = 1.61723e-05
 Nonlinear solver converged at step 3


*** Solving time step 14, time = 0.15 ***
Linear solver converged at step: 8, final residual: 0.000159547  Nonlinear convergence: ||u - u_old|| = 0.00594451
Linear solver converged at step: 17, final residual: 5.39693e-06  Nonlinear convergence: ||u - u_old|| = 0.00369967
Linear solver converged at step: 250, final residual: 4.58011e-08  Nonlinear convergence: ||u - u_old|| = 0.00112625
Linear solver converged at step: 250, final residual: 7.19345e-10  Nonlinear convergence: ||u - u_old|| = 7.43234e-05
 Nonlinear solver converged at step 3

-------------------------------------------------------
| Time:           Wed Jun  6 12:19:06 2007             |
| OS:             Linux                                |
| HostName:       orville                              |
| OS Release:     2.6.21-1.3194.fc7PAE                 |
| OS Version:     #1 SMP Wed May 23 22:27:31 EDT 2007  |
| Machine:        i686                                 |
| Username:       peterson                             |
-------------------------------------------------------
 ------------------------------------------------------------------------------
| Example 13 Performance: Alive time=30.8237, Active time=30.4748              |
 ------------------------------------------------------------------------------
| Event                         nCalls    Total       Avg         Percent of   |
|                                         Time        Time        Active Time  |
|------------------------------------------------------------------------------|
|                                                                              |
| linear solve                  61        30.4748     0.499587    100.00       |
 ------------------------------------------------------------------------------
| Totals:                       61        30.4748                 100.00       |
 ------------------------------------------------------------------------------

 
***************************************************************
* Done Running Example  ./ex13-devel
***************************************************************