# drivencavity-compatible.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.

The script is identical to drivencavity.py except that it uses the Raviart-Thomas discretization providing compatible velocity and pressure spaces resulting in a pointwise divergence-free velocity field.

```12from nutils import mesh, function, solver, export, cli, testing
13from nutils.expression_v2 import Namespace
14import numpy
15import treelog
16
17
18def main(nelems: int, degree: int, reynolds: float):
19    '''
20    Driven cavity benchmark problem using compatible spaces.
21
22    .. arguments::
23
24       nelems [12]
25         Number of elements along edge.
26       degree [2]
27         Polynomial degree for velocity; the pressure space is one degree less.
28       reynolds [1000]
29         Reynolds number, taking the domain size as characteristic length.
30    '''
31
32    verts = numpy.linspace(0, 1, nelems+1)
33    domain, geom = mesh.rectilinear([verts, verts])
34
35    ns = Namespace()
36    ns.δ = function.eye(domain.ndims)
37    ns.Σ = function.ones([domain.ndims])
38    ns.x = geom
39    ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
40    ns.Re = reynolds
41    ns.ubasis = function.vectorize([
42        domain.basis('spline', degree=(degree, degree-1), removedofs=((0, -1), None)),
43        domain.basis('spline', degree=(degree-1, degree), removedofs=(None, (0, -1)))])
44    ns.pbasis = domain.basis('spline', degree=degree-1)
45    ns.u = function.dotarg('u', ns.ubasis)
46    ns.p = function.dotarg('p', ns.pbasis)
47    ns.lm = function.Argument('lm', ())
48    ns.stress_ij = '(∇_j(u_i) + ∇_i(u_j)) / Re - p δ_ij'
49    ns.uwall = function.stack([domain.boundary.indicator('top'), 0])
50    ns.N = 5 * degree * nelems  # nitsche constant based on element size = 1/nelems
51    ns.nitsche_ni = '(N ubasis_ni - (∇_j(ubasis_ni) + ∇_i(ubasis_nj)) n_j) / Re'
52
53    ures = domain.integral('∇_j(ubasis_ni) stress_ij dV' @ ns, degree=2*degree)
54    ures += domain.boundary.integral('(nitsche_ni (u_i - uwall_i) - ubasis_ni stress_ij n_j) dS' @ ns, degree=2*degree)
55    pres = domain.integral('pbasis_n (∇_k(u_k) + lm) dV' @ ns, degree=2*degree)
56    lres = domain.integral('p dV' @ ns, degree=2*degree)
57
58    with treelog.context('stokes'):
59        state0 = solver.solve_linear(['u', 'p', 'lm'], [ures, pres, lres])
60        postprocess(domain, ns, **state0)
61
62    ures += domain.integral('ubasis_ni ∇_j(u_i) u_j dV' @ ns, degree=3*degree)
63    with treelog.context('navierstokes'):
64        state1 = solver.newton(('u', 'p', 'lm'), (ures, pres, lres), arguments=state0).solve(tol=1e-10)
65        postprocess(domain, ns, **state1)
66
67    return state0, state1
```

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.

```74def postprocess(domain, ns, every=.05, spacing=.01, **arguments):
75
76    div = domain.integral('∇_k(u_k)^2 dV' @ ns, degree=1).eval(**arguments)**.5
77    treelog.info('velocity divergence: {:.2e}'.format(div))  # confirm that velocity is pointwise divergence-free
78
79    ns = ns.copy_()  # copy namespace so that we don't modify the calling argument
80    ns.streambasis = domain.basis('std', degree=2)[1:]  # remove first dof to obtain non-singular system
81    ns.stream = function.dotarg('streamdofs', ns.streambasis)  # stream function
82    ns.ε = function.levicivita(2)
83    sqr = domain.integral('Σ_i (u_i - ε_ij ∇_j(stream))^2 dV' @ ns, degree=4)
84    arguments['streamdofs'] = solver.optimize('streamdofs', sqr, arguments=arguments)  # compute streamlines
85
86    bezier = domain.sample('bezier', 9)
87    x, u, p, stream = bezier.eval(['x_i', 'sqrt(u_i u_i)', 'p', 'stream'] @ ns, **arguments)
88    with export.mplfigure('flow.png') as fig:  # plot velocity as field, pressure as contours, streamlines as dashed
89        ax = fig.add_axes([.1, .1, .8, .8], yticks=[], aspect='equal')
90        import matplotlib.collections
92        ax.tricontour(x[:, 0], x[:, 1], bezier.tri, stream, 16, colors='k', linestyles='dotted', linewidths=.5, zorder=9)
93        caxu = fig.add_axes([.1, .1, .03, .8], title='velocity')
94        imu = ax.tripcolor(x[:, 0], x[:, 1], bezier.tri, u, shading='gouraud', cmap='jet')
95        fig.colorbar(imu, cax=caxu)
96        caxu.yaxis.set_ticks_position('left')
97        caxp = fig.add_axes([.87, .1, .03, .8], title='pressure')
98        imp = ax.tricontour(x[:, 0], x[:, 1], bezier.tri, p, 16, cmap='gray', linestyles='solid')
99        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-compatible.py``` (view log).

```107if __name__ == '__main__':
108    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.

```117class test(testing.TestCase):
118
119    def test_p2(self):
120        state0, state1 = main(nelems=3, reynolds=100, degree=2)
121        with self.subTest('stokes-velocity'):
122            self.assertAlmostEqual64(state0['u'], '''
123                eNrzu9Bt8OuUndkD/eTTSqezzP2g/E3698/ZmZlf2GjSaHJS3/90/Wm/C4qGh04CAErOF+g=''')
124        with self.subTest('stokes-pressure'):
125            self.assertAlmostEqual64(state0['p'], '''
127        with self.subTest('stokes-multiplier'):
128            self.assertAlmostEqual(state0['lm'], 0)
129        with self.subTest('navier-stokes-velocity'):
130            self.assertAlmostEqual64(state1['u'], '''
132        with self.subTest('navier-stokes-pressure'):
133            self.assertAlmostEqual64(state1['p'], '''
134                eNpbqpasrWa46TSTfoT+krO/de/pntA3Pnf2zFNtn2u2hglmAOKVDlE=''')
135        with self.subTest('navier-stokes-multiplier'):
136            self.assertAlmostEqual(state1['lm'], 0)
137
138    def test_p3(self):
139        state0, state1 = main(nelems=3, reynolds=100, degree=3)
140        with self.subTest('stokes-velocity'):
141            self.assertAlmostEqual64(state0['u'], '''
142                eNqbo3f7/AeDb2dmGGtd7r+odk7icoQJgjUHLpty8b/BvrMzjF/rll9qMD5kwAAFopc6dRvO2J2fo8d4
143                3sko4wwAjL4lyw==''', atol=1e-14)
144        with self.subTest('stokes-pressure'):
145            self.assertAlmostEqual64(state0['p'], '''
146                eNpbrN+us+Z8qWHMyfcXrly013l1ZocRAxwI6uvoHbwsZuxxNvZC5eUQg+5zS8wAElAT9w==''')
147        with self.subTest('stokes-multiplier'):
148            self.assertAlmostEqual(state0['lm'], 0)
149        with self.subTest('navier-stokes-velocity'):
150            self.assertAlmostEqual64(state1['u'], '''
151                eNoBUACv/14yGcxyNPbJYTahLj/LSDE7yy43SM9WMsXJoDR+N3Iw8s1hM5zJODeizcE0X8phNrQwUDOO
152                NbMzJi+ty4s1oDFqzxIzysjWzXIwFM3tNMjIpDMlJw==''')
153        with self.subTest('navier-stokes-pressure'):
154            self.assertAlmostEqual64(state1['p'], '''
155                eNoBMgDN/yoxvDHizaEy88kDLgoviCt4znIzMy7jL7nTRMzqzhAwBTA/LE0w6czY2GQwoNEYMHE3NDUW
156                DQ==''')
157        with self.subTest('navier-stokes-multiplier'):
158            self.assertAlmostEqual(state1['lm'], 0)
```