drivencavity.py¶
In this script we solve the lid driven cavity problem for stationary Stokes and Navier-Stokes flow. That is, a unit square domain, with no-slip left, bottom and right boundaries and a top boundary that is moving at unit velocity in positive x-direction.
8from nutils import mesh, function, solver, export, cli, testing
9import numpy, treelog
10
11def main(nelems:int, etype:str, degree:int, reynolds:float):
12 '''
13 Driven cavity benchmark problem.
14
15 .. arguments::
16
17 nelems [12]
18 Number of elements along edge.
19 etype [square]
20 Element type (square/triangle/mixed).
21 degree [2]
22 Polynomial degree for velocity; the pressure space is one degree less.
23 reynolds [1000]
24 Reynolds number, taking the domain size as characteristic length.
25 '''
26
27 domain, geom = mesh.unitsquare(nelems, etype)
28
29 ns = function.Namespace()
30 ns.Re = reynolds
31 ns.x = geom
32 ns.ubasis, ns.pbasis = function.chain([
33 domain.basis('std', degree=degree).vector(2),
34 domain.basis('std', degree=degree-1),
35 ])
36 ns.u_i = 'ubasis_ni ?lhs_n'
37 ns.p = 'pbasis_n ?lhs_n'
38 ns.stress_ij = '(u_i,j + u_j,i) / Re - p δ_ij'
39
40 sqr = domain.boundary.integral('u_k u_k d:x' @ ns, degree=degree*2)
41 wallcons = solver.optimize('lhs', sqr, droptol=1e-15)
42
43 sqr = domain.boundary['top'].integral('(u_0 - 1)^2 d:x' @ ns, degree=degree*2)
44 lidcons = solver.optimize('lhs', sqr, droptol=1e-15)
45
46 cons = numpy.choose(numpy.isnan(lidcons), [lidcons, wallcons])
47 cons[-1] = 0 # pressure point constraint
48
49 res = domain.integral('(ubasis_ni,j stress_ij + pbasis_n u_k,k) d:x' @ ns, degree=degree*2)
50 with treelog.context('stokes'):
51 lhs0 = solver.solve_linear('lhs', res, constrain=cons)
52 postprocess(domain, ns, lhs=lhs0)
53
54 res += domain.integral('.5 (ubasis_ni u_i,j - ubasis_ni,j u_i) u_j d:x' @ ns, degree=degree*3)
55 with treelog.context('navierstokes'):
56 lhs1 = solver.newton('lhs', res, lhs0=lhs0, constrain=cons).solve(tol=1e-10)
57 postprocess(domain, ns, lhs=lhs1)
58
59 return lhs0, lhs1
Postprocessing in this script is separated so that it can be reused for the results of Stokes and Navier-Stokes, and because of the extra steps required for establishing streamlines.
65def postprocess(domain, ns, every=.05, spacing=.01, **arguments):
66
67 ns = ns.copy_() # copy namespace so that we don't modify the calling argument
68 ns.streambasis = domain.basis('std', degree=2)[1:] # remove first dof to obtain non-singular system
69 ns.stream = 'streambasis_n ?streamdofs_n' # stream function
70 sqr = domain.integral('((u_0 - stream_,1)^2 + (u_1 + stream_,0)^2) d:x' @ ns, degree=4)
71 arguments['streamdofs'] = solver.optimize('streamdofs', sqr, arguments=arguments) # compute streamlines
72
73 bezier = domain.sample('bezier', 9)
74 x, u, p, stream = bezier.eval(['x_i', 'sqrt(u_k u_k)', 'p', 'stream'] @ ns, **arguments)
75 with export.mplfigure('flow.png') as fig: # plot velocity as field, pressure as contours, streamlines as dashed
76 ax = fig.add_axes([.1,.1,.8,.8], yticks=[], aspect='equal')
77 import matplotlib.collections
78 ax.add_collection(matplotlib.collections.LineCollection(x[bezier.hull], colors='w', linewidths=.5, alpha=.2))
79 ax.tricontour(x[:,0], x[:,1], bezier.tri, stream, 16, colors='k', linestyles='dotted', linewidths=.5, zorder=9)
80 caxu = fig.add_axes([.1,.1,.03,.8], title='velocity')
81 imu = ax.tripcolor(x[:,0], x[:,1], bezier.tri, u, shading='gouraud', cmap='jet')
82 fig.colorbar(imu, cax=caxu)
83 caxu.yaxis.set_ticks_position('left')
84 caxp = fig.add_axes([.87,.1,.03,.8], title='pressure')
85 imp = ax.tricontour(x[:,0], x[:,1], bezier.tri, p, 16, cmap='gray', linestyles='solid')
86 fig.colorbar(imp, cax=caxp)
If the script is executed (as opposed to imported), nutils.cli.run()
calls the main function with arguments provided from the command line. To
keep with the default arguments simply run python3 drivencavity.py
(view log).
92if __name__ == '__main__':
93 cli.run(main)
Once a simulation is developed and tested, it is good practice to save a few
strategic return values for regression testing. The nutils.testing
module, which builds on the standard unittest
framework, facilitates
this by providing nutils.testing.TestCase.assertAlmostEqual64()
for the
embedding of desired results as compressed base64 data.
101class test(testing.TestCase):
102
103 @testing.requires('matplotlib')
104 def test_square(self):
105 lhs0, lhs1 = main(nelems=3, etype='square', reynolds=100, degree=3)
106 with self.subTest('stokes'): self.assertAlmostEqual64(lhs0, '''
107 eNqNkE8og3EYx99JLZKDJTktOfjT9nr3+3kjuTi4zcUB5cpcxcGkOa3kQC4OGyUXaRfESiI7OdjzfX/v
108 3nd5OUiWhpbUWlsr/1Za48Tzrec5PJ/69v1K0t8z3PN9F11TZtQsdPmS9WzaauWW/sH9iaNuRe9TbayO
109 lblHkXFFtbzSqU2wNJq4E588JZZ5xNPGvepLoszta+ftGbiVAJY8/RhhaSqymJFkZ+yQhVXLXeYKWOkY
110 RVbOk6yc0J2yQ2MeSX5TgnxVuVcnf3CzV6OoT+TJECfUInZoV5PkahHkM+Je3TAqvgNWBqYIYF7rRwRp
111 siNmuHDGhhBWO4xKjkYzqtWKTm0TaTyTEzZKiTmKeG7IqzrkSi8hV9Ss0X3JLKatW7L0KvInvHFFv7i0
112 sRzK3H96TtErPVGKArQdD8Wv6YEMOqUFGqM9uqNM6WNRjLbomHK/VIv3UgoHZIyjFw2oRjM41nGNQdhR
113 JFtpr+HAlKQvrHfV6g==''')
114 with self.subTest('navier-stokes'): self.assertAlmostEqual64(lhs1, '''
115 eNpjYCAMgswh9ErtA5c9ruTp/7xiZbhX28jo8MVbRn1XLIxVLhcbBxuIGsDUHbhgoxNx3tnA9vwcQ4fz
116 PkbC59mNP2rONOIxnGiUaOx9GaYu5PxOjfJzB/Xdz5kbWp9jNYo9t81ont4so1OGXUbNJtX6MHVd515p
117 XT4rqT//7EWDprPzDR+eTTY6pBdltNhoq9ELE+3zMHWe58rVl53lv2R1Vvbq+zOTdCLPmhpONkgwlDMO
118 NawyvWAIU/f7nMRN03PyF3rOSp6XOqt3bv+ZA+cirnSclzf2vxhvGgo3T/pC5xXW80Ln751ddzb1rOpZ
119 wzOXTp89l3iu1vic0WrTImOYugCd8uvrriy/aHf1w9nkizPPvDl/7rTe+cRTBmf3nVQx0T4DU0dMOK8/
120 vfX0xtOrT3ef9jitfXrN6Q1A9oLTk0/POL339LrTW4AiC05POt0F5G85vR0oMut0z+mK04tP7wfK7zi9
121 /nTf6aWnt50+enrP6ROnL57+cXrfacYze4G8rUD843SpLgMDAGbLx+o=''')
122
123 @testing.requires('matplotlib')
124 def test_mixed(self):
125 lhs0, lhs1 = main(nelems=3, etype='mixed', reynolds=100, degree=2)
126 with self.subTest('stokes'): self.assertAlmostEqual64(lhs0, '''
127 eNpjYEAFEy++uHzs/EYjEFv73Hm9T2eVDGFyMWdfGJ/TVTGxMvQyqjxbAlYTZM7AsNWIxwhEG5hs0wTR
128 1Wd5DDTOHrx05OzmczC9Cud8riSccbrqYJR2dfMZ07M+hkznLpuongepB+Eyk/TzIHVbL4QbrLqw9Qyy
129 m+aee31awWALXEzP+NIZkB6Y/TFns88mn008u/mM+dnYM0Fn28/OOnv27K6zK8+Wn10FdAEAKXRKgw==''')
130 with self.subTest('navier-stokes'): self.assertAlmostEqual64(lhs1, '''
131 eNpjYEAFvRcfXLa8eNsQxP5xLkBv17ktBjC5iHOfjbr1XxotNbAy0j4nbQQSCzJnYLA0lNIH0XLGi3VB
132 dNK5HoNjZ2frXj+77BxM775zfbq6Z+cYxBpd1vc543jursGSC0omz86D1IPwBJOzYDszL5oaLbgQcQbZ
133 TTZnec9svix4FsZXNbl9GqQHZj/rGb4zPGfYz5w5/fr03tOMZzjOKJz5d5rzjMmZL6cvAk0CAOBkR7A=''')