# 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
79    ax.tricontour(x[:,0], x[:,1], bezier.tri, stream, 16, colors='k', linestyles='dotted', linewidths=.5, zorder=9)
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')
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=''')
```