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=''')