NSASM
Sparse Navier-Stokes Jacobian Assembly

NSASM is a library of C routines which carry 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.

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 is 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 ) ) );
        end
      

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:

        N3
         |\
         | \
         |  \
        N5   N4
         |    \
         |     \
         |      \
        N1--N6--N2
      

Usage:

        [ K, L ] = nsasm ( p, t, np0, e, u, mu );
      

real P[3,NP]

the coordinates of the nodes;

integer T[6,NT]

the indices of the nodes that make up the elements; the node indices are 1-based, not 0-based; the ordering of the 6 entries is vertex 1, vertex 2, vertex 3, then the node between vertex 2 and vertex 3, then the node between vertex 3 and vertex 1, and then the node between vertex 1 and vertex 2. (This is a somewhat peculiar ordering.)

integer NP0

the number of pressure nodes, that is, nodes which are not the vertices of a triangular element.

real E[3,NE]

the contstraints. Each constraint is defined by listing a node index, a variable type (0=horizonal velocity, 1 = vertical velocity, 2 = pressure), and a value to be assigned that variable.

real U[2*NP+NP0+NE]

the current estimate of the solution.

real MU

the kinematic viscosity.

sparse K

a sparse matrix containing the Jacobian matrix for the Newton iteration.

real L[2*NP+NP0+NE]

the residual vector.

Related Data and Programs:

BUMP is an executable FORTRAN90 program which solves a fluid flow problem in a channel including a bump which obstructs and redirects the flow.

DIRECTION_ARROWS_GRID is a MATLAB program which reads files of node and velocity data, and, using interpolation, creates a velocity direction plot with arrows place on a uniform grid of the user's specification.

DISTMESH is a MATLAB program which computes an unstructured mesh.

FEM1D, is an executable FORTRAN77 program which applies the finite element method, with piecewise linear basis functions, to a two point boundary value problem;

FREE_FEM_NAVIER_STOKES is a MATLAB program for solving the steady incompressible Navier Stokes equations on an arbitrary triangulated region, using the finite element method.

MATLAB_C contains some demonstrations of how a MATLAB procedure can invoke an underlying C function.

SPARSE describes MATLAB's sparse function, which can be used to define or access a sparse matrix.

SUPERLU is a library of C routines which can carry out a fast direct solution of large sparse linear systems.

VECTOR_PLOT is an executable FORTRAN90 program that can plot the velocity field and the velocity direction field.

VECTOR_STREAM_GRID is a MATLAB program which reads node and vector data from a file, computes an interpolatory function, evaluates on a uniform grid of points specified by the user, and displays a streamline plot of the vector field.

VELOCITY_ARROWS. is a MATLAB program that can plot the velocity field.

VELOCITY_ARROWS_GRID is a MATLAB program which reads files of node and velocity data, and, using interpolation, creates a vector plot with arrows place on a uniform grid of the user's specification.

Author:

Per-Olof Persson, whose web site is: http://math.mit.edu/~persson/.

Reference:

1. Klaus-Juergen Bathe,
   Finite Element Procedures,
   Prentice Hall, 1996.
2. James Demmel, John Gilbert, Xiaoye Li,
   SuperLU User's Guide,
   September, 1999.
3. Per-Olof Persson,
   Implementation of Finite Element-Based Navier-Stokes Solver,
   April 2002,
   persson_nsasm.pdf
4. Michael Schaefer, Stefan Turek,
   Benchmark Computations of Laminar Flow Around a Cylinder,
   Notes on Numerical Fluid Mechanics,
   Volume 52, 1996, pages 547-566.

Source Code:

• nsasm.c, the source code.
• nsasm.h, the include file.
• nsasm.csh, commands to compile the source code as a simple C file, NOT using Matlab's MEX compiler. This is simply to check for syntax errors.

Examples and Tests:

NSASM_INTERFACE is a MATLAB script that reads a set of point, element, and constraint files, and calls NSASM to determine the stiffness matrix and right hand side.

• dsp_print_some.m, a MATLAB routine that prints part of a sparse triplet matrix.
• nsasm_interface.m, a MATLAB routine that calls NSASM.
• s_len_trim.m, a MATLAB routine that returns the length of a string to the last nonblank.
• timestamp.m, a MATLAB routine that prints the current YMDHMS date as a timestamp.

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

• nsasm_big_test.m, calls NSASM_INTERFACE for the "big" data.
• nsasm_big_test.out, the printed output.
• big_constant.txt, records the constants for the problem, NP, NT, NE, NP0, MU.
• big_e.txt, the constraint file (3*1287 entries).
• big_p.txt, the point coordinate file (2*2049 entries).
• big_t.txt, the triangle file (6*960 entries).

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

• nsasm_small_test.m, calls NSASM_INTERFACE for the "small" data.
• nsasm_small_test.out, the printed output.
• small_e.txt, the constraint file (3*33 entries).
• small_p.txt, the point coordinate file (2*25 entries).
• small_t.txt, the triangle file (6*8 entries).

List of Routines:

• mexFunction is the interface between MATLAB and the C code.
• ASSEMBLE assembles the local stiffness and residual into global arrays.
• ASSEMBLE_CONSTR assembles the constraints.
• INIT_SHAPE evaluates the shape functions at the quadrature points.
• IVEC_HEAP_D reorders an array of integers into a descending heap.
• IVEC_SORT_HEAP_A ascending sorts an array of integers using heap sort.
• LOCALKL assembles the local stiffness matrix and residual.
• SPARSE_CREATE creates a sparse matrix.
• SPARSE_SET increments an entry of the sparse matrix.

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