-
Notifications
You must be signed in to change notification settings - Fork 59
/
Copy pathnsasm.html
314 lines (273 loc) · 8.53 KB
/
nsasm.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
<html>
<head>
<title>
NSASM - Sparse Navier-Stokes Jacobian Assembly
</title>
</head>
<body bgcolor="#EEEEEE" link="#CC0000" alink="#FF3300" vlink="#000055">
<h1 align = "center">
NSASM <br> Sparse Navier-Stokes Jacobian Assembly
</h1>
<hr>
<p>
<b>NSASM</b>
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.
</p>
<p>
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.
</p>
<p>
The sparse matrix format used is known as the <i>compressed
column format</i>, in which the nonzero entries and their row indices
are stored in order by columns.
</p>
<p>
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.
</p>
<p>
The software was originally designed to be called,
via MATLAB's <b>mex</b> facility,
with the linear system solved by <b>SUPER_LU</b>.
Thus, the Newton iteration could have the following form:
<pre>
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
</pre>
</p>
<h3 align = "center">
Indexing Conventions
</h3>
<p>
The software makes some particular assumptions about the
numbering of nodes, variables, and local element nodes
</p>
<p>
<i>(Indexing of nodes)</i> 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.
</p>
<p>
<i>(Indexing of global variables)</i> 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.
</p>
<p>
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.
</p>
<p>
<i>(Numerical codes for U,V,P)</i> In the E vector used to define
constraints, the numeric codes 0, 1 and 2 are used for U, V and
P variables respectively.
</p>
<p>
<i>(Local numbering of nodes)</i> When listing the 6 nodes
that make up a triangle, the ordering should be as follows:
<pre>
N3
|\
| \
| \
N5 N4
| \
| \
| \
N1--N6--N2
</pre>
</p>
<h3 align = "center">
Licensing:
</h3>
<p>
The computer code and data files described and made available on this
web page are distributed under
<a href = "../../txt/gnu_lgpl.txt">the GNU LGPL license.</a>
</p>
<h3 align = "center">
Languages:
</h3>
<p>
<b>NSASM</b> is available in
<a href = "../../c_src/nsasm/nsasm.html"> a C version</a> and
<a href = "../../f_src/nsasm/nsasm.html"> a FORTRAN90 version</a>,
</p>
<h3 align = "center">
Related Data and Programs:
</h3>
<p>
<a href = "../../f_src/fem2d_navier_stokes/fem2d_navier_stokes.html">
FEM2D_NAVIER_STOKES</a>,
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.
</p>
<h3 align = "center">
Author:
</h3>
<p>
Original C version by Per-Olof Persson;<br>
FORTRAN90 version by John Burkardt.
</p>
<h3 align = "center">
Reference:
</h3>
<p>
<ol>
<li>
Per-Olof Persson,<br>
Implementation of Finite Element-Based Navier-Stokes Solver,<br>
April 2002.
</li>
</ol>
</p>
<h3 align = "center">
Source Code:
</h3>
<p>
<ul>
<li>
<a href = "nsasm.f90">nsasm.f90</a>, the source code.
</li>
<li>
<a href = "nsasm.sh">nsasm.sh</a>,
BASH commands to compile the source code.
</li>
</ul>
</p>
<h3 align = "center">
Examples and Tests:
</h3>
<p>
The <b>BIG</b> test defines a relatively large sparse matrix,
with 2049 nodes, 960 elements, 1287 constraints, NP0 = 545,
MU = 500.
<ul>
<li>
<a href = "big_constant.txt">big_constant.txt</a>,
records the constants for the problem, NP, NT, NE, NP0, MU.
</li>
<li>
<a href = "big_e.txt">big_e.txt</a>,
the constraint file (3*1287 entries).
</li>
<li>
<a href = "big_nodes.txt">big_nodes.txt</a>,
the point coordinate file, 2049 nodes.
</li>
<li>
<a href = "big_nodes.png">big_nodes.png</a>,
a PNG image of the nodes.
</li>
<li>
<a href = "big_elements.txt">big_elements.txt</a>,
the element file, 960 quadratic (6 node) triangles.
</li>
<li>
<a href = "big_elements.png">big_elements.png</a>,
a PNG image of the elements.
</li>
</ul>
</p>
<p>
The <b>SMALL</b> test defines a small sparse matrix,
with 25 nodes, 8 elements, 33 constraints, NP0 = 9,
MU = 100.
<ul>
<li>
<a href = "small_e.txt">small_e.txt</a>,
the constraint file (3*33 entries).
</li>
<li>
<a href = "small_nodes.txt">small_nodes.txt</a>,
the point coordinate file, 25 nodes.
</li>
<li>
<a href = "small_nodes.png">small_nodes.png</a>,
a PNG image of the nodes.
</li>
<li>
<a href = "small_elements.txt">small_elements.txt</a>,
the element file, 8 quadratic (6 node) triangles.
</li>
<li>
<a href = "small_elements.png">small_elements.png</a>,
a PNG image of the elements.
</li>
</ul>
</p>
<h3 align = "center">
List of Routines:
</h3>
<p>
<ul>
<li>
<b>ASSEMBLE</b> assembles the local stiffness and residual into global arrays.
</li>
<li>
<b>ASSEMBLE_CONSTR</b> assembles the constraints.
</li>
<li>
<b>INIT_SHAPE</b> evaluates the shape functions at the quadrature points.
</li>
<li>
<b>I4VEC_HEAP_D</b> reorders an I4VEC into an descending heap.
</li>
<li>
<b>I4VEC_SORT_HEAP_A</b> ascending sorts an I4VEC using heap sort.
</li>
<li>
<b>LOCALKL</b> assembles the local stiffness matrix and residual.
</li>
<li>
<b>QUAD_RULE</b> returns the points and weights of a quadrature rule.
</li>
<li>
<b>R8SP_PRINT_SOME</b> prints some of an R8SP matrix.
</li>
<li>
<b>SPARSE_SET</b> increments an entry of the sparse matrix.
</li>
<li>
<b>TIMESTAMP</b> prints the current YMDHMS date as a time stamp.
</li>
</ul>
</p>
<p>
You can go up one level to <a href = "../f_src.html">
the FORTRAN90 source codes</a>.
</p>
<hr>
<i>
Last revised on 23 January 2011.
</i>
<!-- John Burkardt -->
</body>
<!-- Initial HTML skeleton created by HTMLINDEX. -->
</html>