NSASM
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 ) ) );
        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
      

Licensing:

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

Languages:

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.

Author:

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

Reference:

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

Source Code:

• nsasm.f90, the 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.

• big_constants.txt, records the constants for the problem, NP, NT, NE, NP0, NU.
• big_constraints.txt, the constraint file (3*1287 entries).
• big_nodes.txt, the point coordinate file, 2049 nodes.
• big_nodes.png, a PNG image of the nodes.
• big_elements.txt, the element file, 960 quadratic (6 node) triangles.
• big_elements.png, a PNG image of the elements.

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

• nsasm_test_small.f90, calls nsasm_interface() for the "small" data.
• small_constants.txt, records the constants for the problem, NP, NT, NE, NP0, NU.
• small_constraints.txt, the constraint file (3*33 entries).
• small_nodes.txt, the point coordinate file, 25 nodes.
• small_nodes.png, a PNG image of the nodes.
• small_elements.txt, the element file, 8 quadratic (6 node) triangles.
• small_elements.png, a PNG image of the elements.

List of Routines:

• ASSEMBLE assembles the local stiffness and residual into global arrays.
• ASSEMBLE_CONSTR assembles the constraints.
• CH_CAP capitalizes a single character.
• CH_EQI is a case insensitive comparison of two characters for equality.
• CH_TO_DIGIT returns the integer value of a base 10 digit.
• FILE_COLUMN_COUNT counts the number of columns in the first line of a file.
• FILE_ROW_COUNT counts the number of row records in a file.
• GET_UNIT returns a free FORTRAN unit number.
• I4MAT_DATA_READ reads data from an I4MAT file.
• I4MAT_HEADER_READ reads the header from an I4MAT.
• I4VEC_HEAP_D reorders an I4VEC into an descending heap.
• I4VEC_SORT_HEAP_A ascending sorts an I4VEC using heap sort.
• INIT_SHAPE evaluates the shape functions at the quadrature points.
• LOCALKL assembles the local stiffness matrix and residual.
• NSASM_INTERFACE is an interface for NSASM.
• QUAD_RULE returns the points and weights of a quadrature rule.
• R8MAT_DATA_READ reads data from an R8MAT file.
• R8MAT_HEADER_READ reads the header from an R8MAT file.
• R8SP_PRINT_SOME prints some of an R8SP matrix.
• R8VEC_PRINT_PART prints "part" of an R8VEC.
• S_TO_R8 reads an R8 from a string.
• S_TO_R8VEC reads an R8VEC from a string.
• S_WORD_COUNT counts the number of "words" in a string.
• SPARSE_SET increments an entry of the sparse matrix.
• TIMESTAMP prints the current YMDHMS date as a time stamp.

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