The source file systems_of_equations_ex2.C with comments:

Systems Example 2 - Unsteady Nonlinear Navier-Stokes



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 <sstream>
        #include <math.h>
        
Basic include file needed for the mesh functionality.
        #include "libmesh/libmesh.h"
        #include "libmesh/mesh.h"
        #include "libmesh/mesh_generation.h"
        #include "libmesh/exodusII_io.h"
        #include "libmesh/equation_systems.h"
        #include "libmesh/fe.h"
        #include "libmesh/quadrature_gauss.h"
        #include "libmesh/dof_map.h"
        #include "libmesh/sparse_matrix.h"
        #include "libmesh/numeric_vector.h"
        #include "libmesh/dense_matrix.h"
        #include "libmesh/dense_vector.h"
        #include "libmesh/linear_implicit_system.h"
        #include "libmesh/transient_system.h"
        #include "libmesh/perf_log.h"
        #include "libmesh/boundary_info.h"
        #include "libmesh/utility.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 "libmesh/dense_submatrix.h"
        #include "libmesh/dense_subvector.h"
        
The definition of a geometric element
        #include "libmesh/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, with dimension to be overridden later, distributed across the default MPI communicator.
          Mesh mesh(init.comm());
        
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("Systems Example 2");
        
Get a reference to the Stokes system to use later.
          TransientLinearImplicitSystem&  navier_stokes_system =
                equation_systems.get_system<TransientLinearImplicitSystem>("Navier-Stokes");
        
Now we begin the timestep loop to compute the time-accurate solution of the equations.
          const Real dt = 0.01;
          navier_stokes_system.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;
        
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 step size as a parameter in the EquationSystem.
              navier_stokes_system.time += dt;
        
A pretty update message
              std::cout << "\n\n*** Solving time step " << t_step <<
                           ", time = " << navier_stokes_system.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)
                {
                  std::ostringstream file_name;
        
We write the file in the ExodusII format.
                  file_name << "out_"
                            << std::setw(3)
                            << std::setfill('0')
                            << std::right
                            << t_step + 1
                            << ".e";
        
                  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_equal_to (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<dof_id_type> dof_indices;
          std::vector<dof_id_type> dof_indices_u;
          std::vector<dof_id_type> dof_indices_v;
          std::vector<dof_id_type> 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 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)
                    {
                      AutoPtr<Elem> side (elem->build_side(s));
        
Loop over the nodes on the side.
                      for (unsigned int ns=0; ns<side->n_nodes(); ns++)
                        {
Boundary ids are set internally by build_square(). 0=bottom 1=right 2=top 3=left

Set u = 1 on the top boundary, 0 everywhere else
                          const Real u_value = (mesh.boundary_info->has_boundary_id(elem,s,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 source file systems_of_equations_ex2.C without comments:

 
  
  #include <iostream>
  #include <algorithm>
  #include <sstream>
  #include <math.h>
  
  #include "libmesh/libmesh.h"
  #include "libmesh/mesh.h"
  #include "libmesh/mesh_generation.h"
  #include "libmesh/exodusII_io.h"
  #include "libmesh/equation_systems.h"
  #include "libmesh/fe.h"
  #include "libmesh/quadrature_gauss.h"
  #include "libmesh/dof_map.h"
  #include "libmesh/sparse_matrix.h"
  #include "libmesh/numeric_vector.h"
  #include "libmesh/dense_matrix.h"
  #include "libmesh/dense_vector.h"
  #include "libmesh/linear_implicit_system.h"
  #include "libmesh/transient_system.h"
  #include "libmesh/perf_log.h"
  #include "libmesh/boundary_info.h"
  #include "libmesh/utility.h"
  
  #include "libmesh/dense_submatrix.h"
  #include "libmesh/dense_subvector.h"
  
  #include "libmesh/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(init.comm());
  
    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("Systems Example 2");
  
    TransientLinearImplicitSystem&  navier_stokes_system =
          equation_systems.get_system<TransientLinearImplicitSystem>("Navier-Stokes");
  
    const Real dt = 0.01;
    navier_stokes_system.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;
  
    AutoPtr<NumericVector<Number> >
      last_nonlinear_soln (navier_stokes_system.solution->clone());
  
    for (unsigned int t_step=0; t_step<n_timesteps; ++t_step)
      {
        navier_stokes_system.time += dt;
  
        std::cout << "\n\n*** Solving time step " << t_step <<
                     ", time = " << navier_stokes_system.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)
          {
            std::ostringstream file_name;
  
            file_name << "out_"
                      << std::setw(3)
                      << std::setfill('0')
                      << std::right
                      << t_step + 1
                      << ".e";
  
            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_equal_to (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<dof_id_type> dof_indices;
    std::vector<dof_id_type> dof_indices_u;
    std::vector<dof_id_type> dof_indices_v;
    std::vector<dof_id_type> 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)
              {
                AutoPtr<Elem> side (elem->build_side(s));
  
                for (unsigned int ns=0; ns<side->n_nodes(); ns++)
                  {
  
                    const Real u_value = (mesh.boundary_info->has_boundary_id(elem,s,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:

make[4]: Entering directory `/net/spark/workspace/roystgnr/libmesh/git/devel/examples/systems_of_equations/systems_of_equations_ex2'
***************************************************************
* Running Example systems_of_equations_ex2:
*  mpirun -np 4 example-devel  -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc
***************************************************************
 
 Mesh Information:
  mesh_dimension()=2
  spatial_dimension()=3
  n_nodes()=1681
    n_local_nodes()=445
  n_elem()=400
    n_local_elem()=100
    n_active_elem()=400
  n_subdomains()=1
  n_partitions()=4
  n_processors()=4
  n_threads()=1
  processor_id()=0

 EquationSystems
  n_systems()=1
   System #0, "Navier-Stokes"
    Type "TransientLinearImplicit"
    Variables={ "u" "v" } "p" 
    Finite Element Types="LAGRANGE", "JACOBI_20_00" "LAGRANGE", "JACOBI_20_00" 
    Infinite Element Mapping="CARTESIAN" "CARTESIAN" 
    Approximation Orders="SECOND", "THIRD" "FIRST", "THIRD" 
    n_dofs()=3803
    n_local_dofs()=1013
    n_constrained_dofs()=0
    n_local_constrained_dofs()=0
    n_vectors()=3
    n_matrices()=1
    DofMap Sparsity
      Average  On-Processor Bandwidth <= 37.4744
      Average Off-Processor Bandwidth <= 2.15514
      Maximum  On-Processor Bandwidth <= 63
      Maximum Off-Processor Bandwidth <= 37
    DofMap Constraints
      Number of DoF Constraints = 0
      Number of Node Constraints = 0



*** Solving time step 0, time = 0.01 ***
Linear solver converged at step: 65, final residual: 0.000298248  Nonlinear convergence: ||u - u_old|| = 282.478
Linear solver converged at step: 25, final residual: 2.45925e-05  Nonlinear convergence: ||u - u_old|| = 0.668719
Linear solver converged at step: 40, final residual: 1.44765e-07  Nonlinear convergence: ||u - u_old|| = 0.0286564
Linear solver converged at step: 83, final residual: 5.42411e-12  Nonlinear convergence: ||u - u_old|| = 0.000141714
 Nonlinear solver converged at step 3


*** Solving time step 1, time = 0.02 ***
Linear solver converged at step: 59, final residual: 0.000259269  Nonlinear convergence: ||u - u_old|| = 30.7567
Linear solver converged at step: 45, final residual: 2.05194e-05  Nonlinear convergence: ||u - u_old|| = 0.0414233
Linear solver converged at step: 39, final residual: 1.19625e-07  Nonlinear convergence: ||u - u_old|| = 0.0328727
Linear solver converged at step: 87, final residual: 4.01951e-12  Nonlinear convergence: ||u - u_old|| = 8.44974e-05
 Nonlinear solver converged at step 3


*** Solving time step 2, time = 0.03 ***
Linear solver converged at step: 53, final residual: 0.000240605  Nonlinear convergence: ||u - u_old|| = 6.11214
Linear solver converged at step: 27, final residual: 1.73995e-05  Nonlinear convergence: ||u - u_old|| = 0.010112
Linear solver converged at step: 51, final residual: 8.64983e-08  Nonlinear convergence: ||u - u_old|| = 0.0208942
Linear solver converged at step: 85, final residual: 2.23364e-12  Nonlinear convergence: ||u - u_old|| = 8.57952e-05
 Nonlinear solver converged at step 3


*** Solving time step 3, time = 0.04 ***
Linear solver converged at step: 47, final residual: 0.000293691  Nonlinear convergence: ||u - u_old|| = 2.15387
Linear solver converged at step: 28, final residual: 2.47274e-05  Nonlinear convergence: ||u - u_old|| = 0.352608
Linear solver converged at step: 42, final residual: 1.32453e-07  Nonlinear convergence: ||u - u_old|| = 0.00767854
Linear solver converged at step: 84, final residual: 5.30133e-12  Nonlinear convergence: ||u - u_old|| = 9.64303e-05
 Nonlinear solver converged at step 3


*** Solving time step 4, time = 0.05 ***
Linear solver converged at step: 28, final residual: 0.000264123  Nonlinear convergence: ||u - u_old|| = 0.986234
Linear solver converged at step: 28, final residual: 1.83163e-05  Nonlinear convergence: ||u - u_old|| = 0.213293
Linear solver converged at step: 57, final residual: 7.41159e-08  Nonlinear convergence: ||u - u_old|| = 0.00400173
Linear solver converged at step: 96, final residual: 1.70644e-12  Nonlinear convergence: ||u - u_old|| = 4.01715e-05
 Nonlinear solver converged at step 3


*** Solving time step 5, time = 0.06 ***
Linear solver converged at step: 26, final residual: 0.000282338  Nonlinear convergence: ||u - u_old|| = 0.511153
Linear solver converged at step: 28, final residual: 2.03238e-05  Nonlinear convergence: ||u - u_old|| = 0.106921
Linear solver converged at step: 51, final residual: 1.26169e-07  Nonlinear convergence: ||u - u_old|| = 0.000901212
 Nonlinear solver converged at step 2


*** Solving time step 6, time = 0.07 ***
Linear solver converged at step: 24, final residual: 0.000267416  Nonlinear convergence: ||u - u_old|| = 0.280512
Linear solver converged at step: 27, final residual: 1.5249e-05  Nonlinear convergence: ||u - u_old|| = 0.0545105
Linear solver converged at step: 41, final residual: 7.08474e-08  Nonlinear convergence: ||u - u_old|| = 0.00276723
Linear solver converged at step: 85, final residual: 1.15412e-12  Nonlinear convergence: ||u - u_old|| = 5.58692e-05
 Nonlinear solver converged at step 3


*** Solving time step 7, time = 0.08 ***
Linear solver converged at step: 23, final residual: 0.000245679  Nonlinear convergence: ||u - u_old|| = 0.160512
Linear solver converged at step: 27, final residual: 1.66287e-05  Nonlinear convergence: ||u - u_old|| = 0.0251193
Linear solver converged at step: 55, final residual: 6.07431e-08  Nonlinear convergence: ||u - u_old|| = 0.00627766
Linear solver converged at step: 86, final residual: 9.15002e-13  Nonlinear convergence: ||u - u_old|| = 1.84448e-05
 Nonlinear solver converged at step 3


*** Solving time step 8, time = 0.09 ***
Linear solver converged at step: 21, final residual: 0.00024834  Nonlinear convergence: ||u - u_old|| = 0.0950392
Linear solver converged at step: 16, final residual: 1.90121e-05  Nonlinear convergence: ||u - u_old|| = 0.00133845
Linear solver converged at step: 32, final residual: 1.11129e-07  Nonlinear convergence: ||u - u_old|| = 0.01663
Linear solver converged at step: 81, final residual: 3.69976e-12  Nonlinear convergence: ||u - u_old|| = 3.90188e-05
 Nonlinear solver converged at step 3


*** Solving time step 9, time = 0.1 ***
Linear solver converged at step: 20, final residual: 0.000228559  Nonlinear convergence: ||u - u_old|| = 0.0579162
Linear solver converged at step: 12, final residual: 1.49823e-05  Nonlinear convergence: ||u - u_old|| = 0.00115626
Linear solver converged at step: 39, final residual: 6.98083e-08  Nonlinear convergence: ||u - u_old|| = 0.00888167
Linear solver converged at step: 88, final residual: 1.27257e-12  Nonlinear convergence: ||u - u_old|| = 4.30319e-05
 Nonlinear solver converged at step 3


*** Solving time step 10, time = 0.11 ***
Linear solver converged at step: 19, final residual: 0.000218903  Nonlinear convergence: ||u - u_old|| = 0.036109
Linear solver converged at step: 12, final residual: 1.09784e-05  Nonlinear convergence: ||u - u_old|| = 0.00104399
Linear solver converged at step: 47, final residual: 3.68532e-08  Nonlinear convergence: ||u - u_old|| = 0.00460834
Linear solver converged at step: 87, final residual: 4.10491e-13  Nonlinear convergence: ||u - u_old|| = 3.42478e-05
 Nonlinear solver converged at step 3


*** Solving time step 11, time = 0.12 ***
Linear solver converged at step: 18, final residual: 0.000267997  Nonlinear convergence: ||u - u_old|| = 0.022915
Linear solver converged at step: 11, final residual: 2.09843e-05  Nonlinear convergence: ||u - u_old|| = 0.0013608
Linear solver converged at step: 50, final residual: 1.26347e-07  Nonlinear convergence: ||u - u_old|| = 0.00214635
Linear solver converged at step: 89, final residual: 4.9936e-12  Nonlinear convergence: ||u - u_old|| = 0.000162734
 Nonlinear solver converged at step 3


*** Solving time step 12, time = 0.13 ***
Linear solver converged at step: 16, final residual: 0.000312613  Nonlinear convergence: ||u - u_old|| = 0.0147624
Linear solver converged at step: 10, final residual: 2.67168e-05  Nonlinear convergence: ||u - u_old|| = 0.0017403
Linear solver converged at step: 43, final residual: 2.14256e-07  Nonlinear convergence: ||u - u_old|| = 0.000984221
 Nonlinear solver converged at step 2


*** Solving time step 13, time = 0.14 ***
Linear solver converged at step: 15, final residual: 0.000296184  Nonlinear convergence: ||u - u_old|| = 0.00981012
Linear solver converged at step: 12, final residual: 2.3183e-05  Nonlinear convergence: ||u - u_old|| = 0.00232355
Linear solver converged at step: 53, final residual: 1.27309e-07  Nonlinear convergence: ||u - u_old|| = 0.000564091
 Nonlinear solver converged at step 2


*** Solving time step 14, time = 0.15 ***
Linear solver converged at step: 13, final residual: 0.000289793  Nonlinear convergence: ||u - u_old|| = 0.00700597
Linear solver converged at step: 14, final residual: 1.88966e-05  Nonlinear convergence: ||u - u_old|| = 0.00302696
Linear solver converged at step: 52, final residual: 1.09943e-07  Nonlinear convergence: ||u - u_old|| = 0.000238694
 Nonlinear solver converged at step 2

 -------------------------------------------------------------------------------------------------------------------
| Processor id:   0                                                                                                 |
| Num Processors: 4                                                                                                 |
| Time:           Fri Apr 19 11:53:16 2013                                                                          |
| OS:             Linux                                                                                             |
| HostName:       spark.ices.utexas.edu                                                                             |
| OS Release:     2.6.32-279.22.1.el6.x86_64                                                                        |
| OS Version:     #1 SMP Tue Feb 5 14:33:39 CST 2013                                                                |
| Machine:        x86_64                                                                                            |
| Username:       roystgnr                                                                                          |
| Configuration:  ../configure  '--enable-everything'                                                               |
|  'METHODS=devel'                                                                                                  |
|  '--prefix=/h2/roystgnr/libmesh-test'                                                                             |
|  'CXX=distcc /usr/bin/g++'                                                                                        |
|  'CC=distcc /usr/bin/gcc'                                                                                         |
|  'FC=distcc /usr/bin/gfortran'                                                                                    |
|  'F77=distcc /usr/bin/gfortran'                                                                                   |
|  'PETSC_DIR=/opt/apps/ossw/libraries/petsc/petsc-3.3-p2'                                                          |
|  'PETSC_ARCH=gcc-system-mkl-gf-10.3.12.361-mpich2-1.4.1p1-cxx-opt'                                                |
|  'SLEPC_DIR=/opt/apps/ossw/libraries/slepc/slepc-3.3-p2-petsc-3.3-p2-cxx-opt'                                     |
|  'TRILINOS_DIR=/opt/apps/ossw/libraries/trilinos/trilinos-10.12.2/sl6/gcc-system/mpich2-1.4.1p1/mkl-gf-10.3.12.361'|
|  'VTK_DIR=/opt/apps/ossw/libraries/vtk/vtk-5.10.0/sl6/gcc-system'                                                 |
|  'HDF5_DIR=/opt/apps/ossw/libraries/hdf5/hdf5-1.8.9/sl6/gcc-system'                                               |
 -------------------------------------------------------------------------------------------------------------------
 -----------------------------------------------------------------------------------------------------------
| Systems Example 2 Performance: Alive time=10.3573, Active time=8.752                                      |
 -----------------------------------------------------------------------------------------------------------
| 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                  56        8.7520      0.156286    8.7520      0.156286    100.00   100.00   |
 -----------------------------------------------------------------------------------------------------------
| Totals:                       56        8.7520                                          100.00            |
 -----------------------------------------------------------------------------------------------------------

 ----------------------------------------------------------------------------------------------------------------
| libMesh Performance: Alive time=10.4609, Active time=10.3765                                                   |
 ----------------------------------------------------------------------------------------------------------------
| 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   |
|----------------------------------------------------------------------------------------------------------------|
|                                                                                                                |
|                                                                                                                |
| DofMap                                                                                                         |
|   add_neighbors_to_send_list()     1         0.0026      0.002613    0.0040      0.003969    0.03     0.04     |
|   build_sparsity()                 1         0.0023      0.002282    0.0051      0.005138    0.02     0.05     |
|   create_dof_constraints()         1         0.0003      0.000264    0.0003      0.000264    0.00     0.00     |
|   distribute_dofs()                1         0.0043      0.004259    0.0103      0.010292    0.04     0.10     |
|   dof_indices()                    27051     0.2784      0.000010    0.2784      0.000010    2.68     2.68     |
|   prepare_send_list()              1         0.0000      0.000035    0.0000      0.000035    0.00     0.00     |
|   reinit()                         1         0.0055      0.005478    0.0055      0.005478    0.05     0.05     |
|                                                                                                                |
| EquationSystems                                                                                                |
|   build_solution_vector()          15        0.0090      0.000599    0.0612      0.004080    0.09     0.59     |
|                                                                                                                |
| ExodusII_IO                                                                                                    |
|   write_nodal_data()               15        1.5292      0.101945    1.5292      0.101945    14.74    14.74    |
|                                                                                                                |
| FE                                                                                                             |
|   compute_shape_functions()        11200     0.0623      0.000006    0.0623      0.000006    0.60     0.60     |
|   init_shape_functions()           112       0.0011      0.000010    0.0011      0.000010    0.01     0.01     |
|                                                                                                                |
| FEMap                                                                                                          |
|   compute_affine_map()             11200     0.0444      0.000004    0.0444      0.000004    0.43     0.43     |
|   init_reference_to_physical_map() 112       0.0030      0.000026    0.0030      0.000026    0.03     0.03     |
|                                                                                                                |
| Mesh                                                                                                           |
|   all_second_order()               1         0.0022      0.002219    0.0022      0.002219    0.02     0.02     |
|   find_neighbors()                 2         0.0024      0.001202    0.0025      0.001264    0.02     0.02     |
|   renumber_nodes_and_elem()        4         0.0004      0.000088    0.0004      0.000088    0.00     0.00     |
|                                                                                                                |
| MeshCommunication                                                                                              |
|   compute_hilbert_indices()        3         0.0053      0.001775    0.0053      0.001775    0.05     0.05     |
|   find_global_indices()            3         0.0008      0.000274    0.0072      0.002394    0.01     0.07     |
|   parallel_sort()                  3         0.0005      0.000179    0.0007      0.000217    0.01     0.01     |
|                                                                                                                |
| MeshOutput                                                                                                     |
|   write_equation_systems()         15        0.0005      0.000033    1.5920      0.106136    0.00     15.34    |
|                                                                                                                |
| MeshTools::Generation                                                                                          |
|   build_cube()                     1         0.0004      0.000417    0.0004      0.000417    0.00     0.00     |
|                                                                                                                |
| MetisPartitioner                                                                                               |
|   partition()                      2         0.0091      0.004545    0.0138      0.006890    0.09     0.13     |
|                                                                                                                |
| Parallel                                                                                                       |
|   allgather()                      11        0.0003      0.000024    0.0003      0.000027    0.00     0.00     |
|   max(bool)                        2         0.0000      0.000007    0.0000      0.000007    0.00     0.00     |
|   max(scalar)                      348       0.0019      0.000005    0.0019      0.000005    0.02     0.02     |
|   max(vector)                      74        0.0005      0.000007    0.0015      0.000021    0.00     0.01     |
|   min(bool)                        412       0.0017      0.000004    0.0017      0.000004    0.02     0.02     |
|   min(scalar)                      341       0.0194      0.000057    0.0194      0.000057    0.19     0.19     |
|   min(vector)                      74        0.0006      0.000008    0.0017      0.000023    0.01     0.02     |
|   probe()                          48        0.0002      0.000003    0.0002      0.000003    0.00     0.00     |
|   receive()                        48        0.0002      0.000004    0.0004      0.000007    0.00     0.00     |
|   send()                           48        0.0001      0.000002    0.0001      0.000002    0.00     0.00     |
|   send_receive()                   54        0.0003      0.000006    0.0009      0.000016    0.00     0.01     |
|   sum()                            65        0.0024      0.000036    0.0192      0.000296    0.02     0.19     |
|                                                                                                                |
| Parallel::Request                                                                                              |
|   wait()                           48        0.0001      0.000001    0.0001      0.000001    0.00     0.00     |
|                                                                                                                |
| Partitioner                                                                                                    |
|   set_node_processor_ids()         2         0.0013      0.000626    0.0020      0.000977    0.01     0.02     |
|   set_parent_processor_ids()       2         0.0003      0.000131    0.0003      0.000131    0.00     0.00     |
|                                                                                                                |
| PetscLinearSolver                                                                                              |
|   solve()                          56        7.2909      0.130194    7.2909      0.130194    70.26    70.26    |
|                                                                                                                |
| System                                                                                                         |
|   assemble()                       56        1.0926      0.019511    1.4521      0.025931    10.53    13.99    |
 ----------------------------------------------------------------------------------------------------------------
| Totals:                            51434     10.3765                                         100.00            |
 ----------------------------------------------------------------------------------------------------------------

 
***************************************************************
* Done Running Example systems_of_equations_ex2:
*  mpirun -np 4 example-devel  -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc
***************************************************************
make[4]: Leaving directory `/net/spark/workspace/roystgnr/libmesh/git/devel/examples/systems_of_equations/systems_of_equations_ex2'

Site Created By: libMesh Developers
Last modified: April 23 2013 04:19:30 UTC

Hosted By:
SourceForge.net Logo