Sparse Navier-Stokes Jacobian Assembly

NSASM is a FORTRAN90 library which carries out the finite element assembly of the jacobian matrix used in a Newton iteration to solve the steady state incompressible Navier Stokes equations in a 2D region, by Per-Olof Persson.

The finite element approximation uses the P2-P1 triangular element, also known as the Taylor-Hood element. The velocities are approximated quadratically, and the pressures linearly.

The sparse matrix format used is known as the compressed column format, in which the nonzero entries and their row indices are stored in order by columns.

The software assumes that a suitable mesh has been set up for the region, with arrays P and T defining that mesh, that an array defining the boundary constraints has been set up, and that an initial estimate of the flow has been made available.

The software was originally designed to be called, via MATLAB's mex facility, with the linear system solved by SUPER_LU. Thus, the Newton iteration could have the following form:

        for ii = 1 : 8
          [ K, L ] = nsasm ( p, t, np0, e, u, mu );
          du = - lusolve ( K, L );
          u = u + du;
          disp ( sprintf ( '%20.10g  %20.10g', norm ( L, inf ), norm ( du, inf ) ) );

Indexing Conventions

The software makes some particular assumptions about the numbering of nodes, variables, and local element nodes

(Indexing of nodes) It is assumed that the mesh was generated by starting with a set of nodes, triangulating them, and then computing the midsides of the triangles and adding these midsides to the list of nodes. In particular, the program therefore assumes that the numbering of the nodes reflects this process, so that all vertex nodes are numbered first, followed by the midside nodes.

(Indexing of global variables) It is assumed that every node has an associated horizontal velocity variable "U", that every node has an associated vertical velocity variable "V", and that only vertex nodes or "pressure nodes" have an associated pressure variable "P". The variables are stored in a single vector, first all the U's, then all the V's, then the P's. Thus variable 1 is the value of U at node 1, variable NP+1 is the value of V at node 1, and variable 2*NP+1 is the value of P at node 1.

Notice that there are even more global variables than you might think, since the boundary conditions and other constraints are handled by adding a Lagrange multiplier for each one. Thus the number of variables is actually 2*NP+NP0+NE.

(Numerical codes for U,V,P) In the E vector used to define constraints, the numeric codes 0, 1 and 2 are used for U, V and P variables respectively.

(Local numbering of nodes) When listing the 6 nodes that make up a triangle, the ordering should be as follows:

         | \
         |  \
        N5   N4
         |    \
         |     \
         |      \


The computer code and data files described and made available on this web page are distributed under the GNU LGPL license.


NSASM is available in a C version and a FORTRAN90 version and a MATLAB version.

Related Data and Programs:

FEM2D_NAVIER_STOKES, a FORTRAN90 program which solves the 2D incompressible Navier Stokes equations in an arbitrary triangulated region. In order to run, it requires user-supplied routines that define problem data.


Original C version by Per-Olof Persson;
FORTRAN90 version by John Burkardt.


  1. Per-Olof Persson,
    Implementation of Finite Element-Based Navier-Stokes Solver,
    April 2002.

Source Code:

Examples and Tests:

The BIG test defines a relatively large sparse matrix, with NP = 2049 nodes, NT = 960 elements, NE = 1287 constraints, NP0 = 545, NU = 500.

The SMALL test defines a small sparse matrix, with NP = 25 nodes, NT = 8 elements, NE = 33 constraints, NP0 = 9, NU = 100.

List of Routines:

You can go up one level to the FORTRAN90 source codes.

Last revised on 23 September 2013.