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 a ExodusII viewer such as ParaView 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 "exodusII_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"

Bring in everything from the libMesh namespace

        using namespace libMesh;

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.

          LibMeshInit init (argc, argv);

Skip this 2D example if libMesh was compiled as 1D-only.

          libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");

Trilinos solver NaNs by default on the zero pressure block. We'll skip this example for now.

          if (libMesh::default_solver_package() == TRILINOS_SOLVERS)
            {
              std::cout << "We skip example 13 when using the Trilinos solvers.\n"
                        << std::endl;
              return 0;
            }

Create a mesh.

          Mesh mesh;

Use the MeshTools::Generation mesh generator to create a uniform 2D grid on the square [-1,1]^2. We instruct the mesh generator to build a mesh of 8x8 \p Quad9 elements in 2D. 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.,
                                               QUAD4);
        
          mesh.all_second_order();

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.push("linear solve");
                  equation_systems.get_system("Navier-Stokes").solve();
                  perf_log.pop("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") =
                    std::min(Utility::pow<2>(final_linear_residual), initial_linear_solver_tol);
        
                } // end nonlinear loop

Write out every nth timestep to file.

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

We write the file in the ExodusII format.

                  file_name << "out_";
                  OSSRealzeroright(file_name,3,0, t_step + 1);
                  file_name << ".exd";
                  
                  ExodusII_IO(mesh).write_equation_systems (file_name.str(),
                                                      equation_systems);
                }
        #endif // #ifdef LIBMESH_HAVE_EXODUS_API
            } // end timestep loop.

All done.

          return 0;
        }

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.

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

Get a constant reference to the mesh object.

          const MeshBase& 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)
                          libmesh_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

If this assembly program were to be used on an adaptive mesh, we would have to apply any hanging node constraint equations

              dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);

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 SparseMatrix::add_matrix() and \p NumericVector::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 "exodusII_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"
  
  using namespace libMesh;
  
  void assemble_stokes (EquationSystems& es,
                        const std::string& system_name);
  
  int main (int argc, char** argv)
  {
    LibMeshInit init (argc, argv);
  
    libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
    
    if (libMesh::default_solver_package() == TRILINOS_SOLVERS)
      {
        std::cout << "We skip example 13 when using the Trilinos solvers.\n"
                  << std::endl;
        return 0;
      }
  
    Mesh mesh;
    
    MeshTools::Generation::build_square (mesh,
                                         20, 20,
                                         0., 1.,
                                         0., 1.,
                                         QUAD4);
  
    mesh.all_second_order();
    
    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.push("linear solve");
            equation_systems.get_system("Navier-Stokes").solve();
            perf_log.pop("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") =
              std::min(Utility::pow<2>(final_linear_residual), initial_linear_solver_tol);
  
          } // end nonlinear loop
        
        const unsigned int write_interval = 1;
        
  #ifdef LIBMESH_HAVE_EXODUS_API
        if ((t_step+1)%write_interval == 0)
          {
            OStringStream file_name;
  
            file_name << "out_";
            OSSRealzeroright(file_name,3,0, t_step + 1);
            file_name << ".exd";
            
            ExodusII_IO(mesh).write_equation_systems (file_name.str(),
                                                equation_systems);
          }
  #endif // #ifdef LIBMESH_HAVE_EXODUS_API
      } // end timestep loop.
    
    return 0;
  }
  
  
  
  
  
  
  void assemble_stokes (EquationSystems& es,
                        const std::string& system_name)
  {
    libmesh_assert (system_name == "Navier-Stokes");
    
    const MeshBase& 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)
                    libmesh_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          
        
        dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
  
        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:

Compiling C++ (in optimized mode) ex13.C...
Linking ex13-opt...
***************************************************************
* Running Example  mpirun -np 2 ./ex13-opt -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc -log_summary
***************************************************************
 
 Mesh Information:
  mesh_dimension()=2
  spatial_dimension()=3
  n_nodes()=1681
    n_local_nodes()=863
  n_elem()=400
    n_local_elem()=200
    n_active_elem()=400
  n_subdomains()=1
  n_processors()=2
  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()=1958
    n_constrained_dofs()=0
    n_vectors()=3



*** Solving time step 0, time = 0.01 ***
Linear solver converged at step: 1, final residual: 0.00952267  Nonlinear convergence: ||u - u_old|| = 282.328
Linear solver converged at step: 0, final residual: 0.00952379  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.00952379  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.00952379  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.00952379  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.00952379  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.00952379  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.00952379  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.00952379  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.00952379  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.00952379  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.00952379  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.00952379  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.00952379  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.00952379  Nonlinear convergence: ||u - u_old|| = 0


*** Solving time step 1, time = 0.02 ***
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0


*** Solving time step 2, time = 0.03 ***
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0


*** Solving time step 3, time = 0.04 ***
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0


*** Solving time step 4, time = 0.05 ***
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0


*** Solving time step 5, time = 0.06 ***
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0


*** Solving time step 6, time = 0.07 ***
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0


*** Solving time step 7, time = 0.08 ***
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0


*** Solving time step 8, time = 0.09 ***
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0


*** Solving time step 9, time = 0.1 ***
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0


*** Solving time step 10, time = 0.11 ***
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0


*** Solving time step 11, time = 0.12 ***
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0


*** Solving time step 12, time = 0.13 ***
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0


*** Solving time step 13, time = 0.14 ***
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0


*** Solving time step 14, time = 0.15 ***
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0
Linear solver converged at step: 0, final residual: 0.0105831  Nonlinear convergence: ||u - u_old|| = 0

-------------------------------------------------------------------
| Processor id:   0                                                |
| Num Processors: 2                                                |
| Time:           Tue Feb 22 12:21:14 2011                         |
| OS:             Linux                                            |
| HostName:       daedalus                                         |
| OS Release:     2.6.32-26-generic                                |
| OS Version:     #46-Ubuntu SMP Tue Oct 26 16:47:18 UTC 2010      |
| Machine:        x86_64                                           |
| Username:       roystgnr                                         |
| Configuration:  ./configure run on Thu Feb 17 16:40:46 CST 2011  |
-------------------------------------------------------------------
 -----------------------------------------------------------------------------------------------------------
| Example 13 Performance: Alive time=92.7848, Active time=92.75                                             |
 -----------------------------------------------------------------------------------------------------------
| Event                         nCalls    Total Time  Avg Time    Total Time  Avg Time    % of Active Time  |
|                                         w/o Sub     w/o Sub     With Sub    With Sub    w/o S    With S   |
|-----------------------------------------------------------------------------------------------------------|
|                                                                                                           |
| linear solve                  225       92.7500     0.412222    92.7500     0.412222    100.00   100.00   |
 -----------------------------------------------------------------------------------------------------------
| Totals:                       225       92.7500                                         100.00            |
 -----------------------------------------------------------------------------------------------------------

************************************************************************************************************************
***             WIDEN YOUR WINDOW TO 120 CHARACTERS.  Use 'enscript -r -fCourier9' to print this document            ***
************************************************************************************************************************

---------------------------------------------- PETSc Performance Summary: ----------------------------------------------

./ex13-opt on a gcc-4.5-l named daedalus with 2 processors, by roystgnr Tue Feb 22 12:21:14 2011
Using Petsc Release Version 3.1.0, Patch 5, Mon Sep 27 11:51:54 CDT 2010

                         Max       Max/Min        Avg      Total 
Time (sec):           9.280e+01      1.00001   9.280e+01
Objects:              1.392e+03      1.00000   1.392e+03
Flops:                4.905e+10      3.14211   3.233e+10  6.465e+10
Flops/sec:            5.285e+08      3.14208   3.483e+08  6.967e+08
MPI Messages:         1.922e+03      1.00000   1.922e+03  3.843e+03
MPI Message Lengths:  1.234e+07      1.01676   6.369e+03  2.448e+07
MPI Reductions:       6.586e+03      1.00000

Flop counting convention: 1 flop = 1 real number operation of type (multiply/divide/add/subtract)
                            e.g., VecAXPY() for real vectors of length N --> 2N flops
                            and VecAXPY() for complex vectors of length N --> 8N flops

Summary of Stages:   ----- Time ------  ----- Flops -----  --- Messages ---  -- Message Lengths --  -- Reductions --
                        Avg     %Total     Avg     %Total   counts   %Total     Avg         %Total   counts   %Total 
 0:      Main Stage: 9.2801e+01 100.0%  6.4654e+10 100.0%  3.843e+03 100.0%  6.369e+03      100.0%  6.116e+03  92.9% 

------------------------------------------------------------------------------------------------------------------------
See the 'Profiling' chapter of the users' manual for details on interpreting output.
Phase summary info:
   Count: number of times phase was executed
   Time and Flops: Max - maximum over all processors
                   Ratio - ratio of maximum to minimum over all processors
   Mess: number of messages sent
   Avg. len: average message length
   Reduct: number of global reductions
   Global: entire computation
   Stage: stages of a computation. Set stages with PetscLogStagePush() and PetscLogStagePop().
      %T - percent time in this phase         %F - percent flops in this phase
      %M - percent messages in this phase     %L - percent message lengths in this phase
      %R - percent reductions in this phase
   Total Mflop/s: 10e-6 * (sum of flops over all processors)/(max time over all processors)
------------------------------------------------------------------------------------------------------------------------
Event                Count      Time (sec)     Flops                             --- Global ---  --- Stage ---   Total
                   Max Ratio  Max     Ratio   Max  Ratio  Mess   Avg len Reduct  %T %F %M %L %R  %T %F %M %L %R Mflop/s
------------------------------------------------------------------------------------------------------------------------

--- Event Stage 0: Main Stage

VecMDot                1 1.0 2.7895e-05 1.9 3.92e+03 1.1 0.0e+00 0.0e+00 1.0e+00  0  0  0  0  0   0  0  0  0  0   273
VecNorm              676 1.0 2.9563e-02 3.5 2.65e+06 1.1 0.0e+00 0.0e+00 6.8e+02  0  0  0  0 10   0  0  0  0 11   174
VecScale             226 1.0 4.4966e-04 1.0 4.43e+05 1.1 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0  1911
VecCopy              467 1.0 1.3275e-03 1.4 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
VecSet               460 1.0 5.8198e-04 1.2 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
VecAXPY              676 1.0 5.4984e-03 1.1 2.65e+06 1.1 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0   935
VecMAXPY               2 1.0 1.2875e-05 1.3 7.83e+03 1.1 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0  1182
VecAssemblyBegin     900 1.0 6.4440e-03 1.0 0.00e+00 0.0 4.5e+02 1.1e+03 2.7e+03  0  0 12  2 41   0  0 12  2 44     0
VecAssemblyEnd       900 1.0 7.1549e-04 1.1 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
VecScatterBegin      676 1.0 3.8245e-03 1.3 0.00e+00 0.0 9.0e+02 1.7e+03 0.0e+00  0  0 23  6  0   0  0 23  6  0     0
VecScatterEnd        676 1.0 3.0619e+0123328.9 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 16  0  0  0  0  16  0  0  0  0     0
VecNormalize         226 1.0 2.6631e-02 4.8 1.33e+06 1.1 0.0e+00 0.0e+00 2.3e+02  0  0  0  0  3   0  0  0  0  4    97
MatMult              226 1.0 3.0646e+011079.8 3.43e+07 1.1 4.5e+02 1.3e+03 0.0e+00 17  0 12  2  0  17  0 12  2  0     2
MatSolve               2 1.0 2.6262e-03 1.1 3.33e+06 1.1 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0  2413
MatLUFactorNum       225 1.0 3.7017e+01 3.1 4.90e+10 3.1 0.0e+00 0.0e+00 0.0e+00 26100  0  0  0  26100  0  0  0  1744
MatILUFactorSym      225 1.0 5.3652e+01 1.1 0.00e+00 0.0 0.0e+00 0.0e+00 2.2e+02 55  0  0  0  3  55  0  0  0  4     0
MatAssemblyBegin     450 1.0 6.1128e-02 2.3 0.00e+00 0.0 6.8e+02 2.2e+04 9.0e+02  0  0 18 60 14   0  0 18 60 15     0
MatAssemblyEnd       450 1.0 7.6824e-02 1.1 0.00e+00 0.0 4.0e+00 3.3e+02 4.6e+02  0  0  0  0  7   0  0  0  0  7     0
MatGetRowIJ          225 1.0 7.1526e-05 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
MatGetOrdering       225 1.0 3.0386e-03 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 4.5e+02  0  0  0  0  7   0  0  0  0  7     0
MatZeroEntries       227 1.0 1.3601e-02 1.3 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
KSPGMRESOrthog         1 1.0 6.1035e-05 1.6 7.83e+03 1.1 0.0e+00 0.0e+00 1.0e+00  0  0  0  0  0   0  0  0  0  0   249
KSPSetup             450 1.0 1.1754e-04 1.1 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
KSPSolve             225 1.0 9.0730e+01 1.0 4.90e+10 3.1 4.5e+02 1.3e+03 1.1e+03 98100 12  2 17  98100 12  2 18   713
PCSetUp              450 1.0 9.0686e+01 1.5 4.90e+10 3.1 0.0e+00 0.0e+00 6.8e+02 81100  0  0 10  81100  0  0 11   712
PCSetUpOnBlocks      225 1.0 9.0685e+01 1.5 4.90e+10 3.1 0.0e+00 0.0e+00 6.8e+02 81100  0  0 10  81100  0  0 11   712
PCApply                2 1.0 2.7211e-03 1.1 3.33e+06 1.1 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0  2329
------------------------------------------------------------------------------------------------------------------------

Memory usage is given in bytes:

Object Type          Creations   Destructions     Memory  Descendants' Mem.
Reports information only for process 0.

--- Event Stage 0: Main Stage

                 Vec    20             20       236912     0
         Vec Scatter   229            229       198772     0
           Index Set   908            908      5762128     0
   IS L to G Mapping     3              3        27252     0
              Matrix   228            228   2249441964     0
       Krylov Solver     2              2        18880     0
      Preconditioner     2              2         1408     0
========================================================================================================================
Average time to get PetscTime(): 0
Average time for MPI_Barrier(): 1.38283e-06
Average time for zero size MPI_Send(): 5.48363e-06
#PETSc Option Table entries:
-ksp_right_pc
-log_summary
-pc_type bjacobi
-sub_pc_factor_levels 4
-sub_pc_factor_zeropivot 0
-sub_pc_type ilu
#End of PETSc Option Table entries
Compiled without FORTRAN kernels
Compiled with full precision matrices (default)
sizeof(short) 2 sizeof(int) 4 sizeof(long) 8 sizeof(void*) 8 sizeof(PetscScalar) 8
Configure run at: Fri Oct 15 13:01:23 2010
Configure options: --with-debugging=false --COPTFLAGS=-O3 --CXXOPTFLAGS=-O3 --FOPTFLAGS=-O3 --with-clanguage=C++ --with-shared=1 --with-mpi-dir=/org/centers/pecos/LIBRARIES/MPICH2/mpich2-1.2.1-gcc-4.5-lucid --with-mumps=true --download-mumps=ifneeded --with-parmetis=true --download-parmetis=ifneeded --with-superlu=true --download-superlu=ifneeded --with-superludir=true --download-superlu_dist=ifneeded --with-blacs=true --download-blacs=ifneeded --with-scalapack=true --download-scalapack=ifneeded --with-hypre=true --download-hypre=ifneeded --with-blas-lib="[/org/centers/pecos/LIBRARIES/MKL/mkl-10.0.3.020-gcc-4.5-lucid/lib/em64t/libmkl_intel_lp64.so,/org/centers/pecos/LIBRARIES/MKL/mkl-10.0.3.020-gcc-4.5-lucid/lib/em64t/libmkl_sequential.so,/org/centers/pecos/LIBRARIES/MKL/mkl-10.0.3.020-gcc-4.5-lucid/lib/em64t/libmkl_core.so]" --with-lapack-lib=/org/centers/pecos/LIBRARIES/MKL/mkl-10.0.3.020-gcc-4.5-lucid/lib/em64t/libmkl_solver_lp64_sequential.a
-----------------------------------------
Libraries compiled on Fri Oct 15 13:01:23 CDT 2010 on atreides 
Machine characteristics: Linux atreides 2.6.32-25-generic #44-Ubuntu SMP Fri Sep 17 20:05:27 UTC 2010 x86_64 GNU/Linux 
Using PETSc directory: /org/centers/pecos/LIBRARIES/PETSC3/petsc-3.1-p5
Using PETSc arch: gcc-4.5-lucid-mpich2-1.2.1-cxx-opt
-----------------------------------------
Using C compiler: /org/centers/pecos/LIBRARIES/MPICH2/mpich2-1.2.1-gcc-4.5-lucid/bin/mpicxx -Wall -Wwrite-strings -Wno-strict-aliasing -O3   -fPIC   
Using Fortran compiler: /org/centers/pecos/LIBRARIES/MPICH2/mpich2-1.2.1-gcc-4.5-lucid/bin/mpif90 -fPIC -Wall -Wno-unused-variable -O3    
-----------------------------------------
Using include paths: -I/org/centers/pecos/LIBRARIES/PETSC3/petsc-3.1-p5/gcc-4.5-lucid-mpich2-1.2.1-cxx-opt/include -I/org/centers/pecos/LIBRARIES/PETSC3/petsc-3.1-p5/include -I/org/centers/pecos/LIBRARIES/PETSC3/petsc-3.1-p5/gcc-4.5-lucid-mpich2-1.2.1-cxx-opt/include -I/org/centers/pecos/LIBRARIES/MPICH2/mpich2-1.2.1-gcc-4.5-lucid/include  
------------------------------------------
Using C linker: /org/centers/pecos/LIBRARIES/MPICH2/mpich2-1.2.1-gcc-4.5-lucid/bin/mpicxx -Wall -Wwrite-strings -Wno-strict-aliasing -O3 
Using Fortran linker: /org/centers/pecos/LIBRARIES/MPICH2/mpich2-1.2.1-gcc-4.5-lucid/bin/mpif90 -fPIC -Wall -Wno-unused-variable -O3  
Using libraries: -Wl,-rpath,/org/centers/pecos/LIBRARIES/PETSC3/petsc-3.1-p5/gcc-4.5-lucid-mpich2-1.2.1-cxx-opt/lib -L/org/centers/pecos/LIBRARIES/PETSC3/petsc-3.1-p5/gcc-4.5-lucid-mpich2-1.2.1-cxx-opt/lib -lpetsc       -lX11 -Wl,-rpath,/org/centers/pecos/LIBRARIES/PETSC3/petsc-3.1-p5/gcc-4.5-lucid-mpich2-1.2.1-cxx-opt/lib -L/org/centers/pecos/LIBRARIES/PETSC3/petsc-3.1-p5/gcc-4.5-lucid-mpich2-1.2.1-cxx-opt/lib -lHYPRE -lsuperlu_dist_2.4 -lcmumps -ldmumps -lsmumps -lzmumps -lmumps_common -lpord -lparmetis -lmetis -lscalapack -lblacs -lsuperlu_4.0 -Wl,-rpath,/org/centers/pecos/LIBRARIES/MKL/mkl-10.0.3.020-gcc-4.5-lucid/lib/em64t -L/org/centers/pecos/LIBRARIES/MKL/mkl-10.0.3.020-gcc-4.5-lucid/lib/em64t -lmkl_solver_lp64_sequential -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lm -Wl,-rpath,/org/centers/pecos/LIBRARIES/MPICH2/mpich2-1.2.1-gcc-4.5-lucid/lib -L/org/centers/pecos/LIBRARIES/MPICH2/mpich2-1.2.1-gcc-4.5-lucid/lib -Wl,-rpath,/org/centers/pecos/LIBRARIES/GCC/gcc-4.5.1-lucid/lib/gcc/x86_64-unknown-linux-gnu/4.5.1 -L/org/centers/pecos/LIBRARIES/GCC/gcc-4.5.1-lucid/lib/gcc/x86_64-unknown-linux-gnu/4.5.1 -Wl,-rpath,/org/centers/pecos/LIBRARIES/GCC/gcc-4.5.1-lucid/lib64 -L/org/centers/pecos/LIBRARIES/GCC/gcc-4.5.1-lucid/lib64 -Wl,-rpath,/org/centers/pecos/LIBRARIES/GCC/gcc-4.5.1-lucid/lib -L/org/centers/pecos/LIBRARIES/GCC/gcc-4.5.1-lucid/lib -ldl -lmpich -lopa -lpthread -lrt -lgcc_s -lmpichf90 -lgfortran -lm -lm -lmpichcxx -lstdc++ -ldl -lmpich -lopa -lpthread -lrt -lgcc_s -ldl  
------------------------------------------
 
***************************************************************
* Done Running Example  mpirun -np 2 ./ex13-opt -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc -log_summary
***************************************************************