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
9from nutils.expression_v2 import Namespace
10import numpy
11import treelog
12
13
14def main(nelems: int, etype: str, degree: int, reynolds: float):
15    '''
16    Driven cavity benchmark problem.
17
18    .. arguments::
19
20       nelems [12]
21         Number of elements along edge.
22       etype [square]
23         Element type (square/triangle/mixed).
24       degree [2]
25         Polynomial degree for velocity; the pressure space is one degree less.
26       reynolds [1000]
27         Reynolds number, taking the domain size as characteristic length.
28    '''
29
30    domain, geom = mesh.unitsquare(nelems, etype)
31
32    ns = Namespace()
33    ns.δ = function.eye(domain.ndims)
34    ns.Σ = function.ones([domain.ndims])
35    ns.Re = reynolds
36    ns.x = geom
37    ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
38    ns.ubasis = domain.basis('std', degree=degree).vector(domain.ndims)
39    ns.pbasis = domain.basis('std', degree=degree-1)
40    ns.u = function.dotarg('u', ns.ubasis)
41    ns.p = function.dotarg('p', ns.pbasis)
42    ns.stress_ij = '(∇_j(u_i) + ∇_i(u_j)) / Re - p δ_ij'
43
44    usqr = domain.boundary.integral('u_k u_k dS' @ ns, degree=degree*2)
45    wallcons = solver.optimize('u', usqr, droptol=1e-15)
46
47    usqr = domain.boundary['top'].integral('(u_0 - 1)^2 dS' @ ns, degree=degree*2)
48    lidcons = solver.optimize('u', usqr, droptol=1e-15)
49
50    ucons = numpy.choose(numpy.isnan(lidcons), [lidcons, wallcons])
51    pcons = numpy.zeros(len(ns.pbasis), dtype=bool)
52    pcons[-1] = True  # constrain pressure to zero in a point
53    cons = dict(u=ucons, p=pcons)
54
55    ures = domain.integral('∇_j(ubasis_ni) stress_ij dV' @ ns, degree=degree*2)
56    pres = domain.integral('pbasis_n ∇_k(u_k) dV' @ ns, degree=degree*2)
57    with treelog.context('stokes'):
58        state0 = solver.solve_linear(('u', 'p'), (ures, pres), constrain=cons)
59        postprocess(domain, ns, **state0)
60
61    ures += domain.integral('.5 (ubasis_ni ∇_j(u_i) - ∇_j(ubasis_ni) u_i) u_j dV' @ ns, degree=degree*3)
62    with treelog.context('navierstokes'):
63        state1 = solver.newton(('u', 'p'), (ures, pres), arguments=state0, constrain=cons).solve(tol=1e-10)
64        postprocess(domain, ns, **state1)
65
66    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.

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

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

```112class test(testing.TestCase):
113
114    def test_square(self):
115        state0, state1 = main(nelems=3, etype='square', reynolds=100, degree=3)
116        with self.subTest('stokes-velocity'):
117            self.assertAlmostEqual64(state0['u'], '''
118                eNpjYCAMgswhdJ1O+uWtl7/rp13hMyq4rmx8/cI/44qLW0wMLliZMhrxGMHUPT//WmfruW8GWudSjJ6d
119                FTeWP/vf+NH5TuNVhurGPqZvL8LUbTi3X+P1WV2D2rPthjZnw4yenflpdODSFaO9RpuNZple14Wp+362
121                H5rOvoSw1/H667OXz9eerTxnc3bV2Wdn2M8euKRzdq+R79lZppqXEP4Qvbz1HNd5rXNzzj47+/KM/FnG
122                M4/Ol59ZZXjzjI+psB4iXGbqbL3MeSHtyqezBdfvnrl+gelMxUWf0wYXjp1iNPpyFqaOmHAGAFFhkvE=''')
123        with self.subTest('stokes-pressure'):
124            self.assertAlmostEqual64(state0['p'], '''
125                eNp7dOb9mRdnHp2pPbPw9MzTN848OXPpzJ4z1Wcizqw/c//Ma6DM9TMHzsw/s+PMFxTIdfbvGfazwmf1
126                zkaftTgrdJblrORZ47NTz94463qW/ezPM4xAcsLZjZcZGAAX7kL6''')
127        with self.subTest('navier-stokes-velocity'):
128            self.assertAlmostEqual64(state1['u'], '''
129                eNpjYCAMgswh9ErtA5c9ruTp/7xiZbhX28jo8MVbRn1XLIxVLhcbBxuIGsDUHbhgoxNx3tnA9vwcQ4fz
130                PkbC59mNP2rONOIxnGiUaOx9GaYu5PxOjfJzB/Xdz5kbWp9jNYo9t81ont4so1OGXUbNJtX6MHVd515p
131                XT4rqT//7EWDprPzDR+eTTY6pBdltNhoq9ELE+3zMHWe58rVl53lv2R1Vvbq+zOTdCLPmhpONkgwlDMO
132                NawyvWAIU/f7nMRN03PyF3rOSp6XOqt3bv+ZA+cirnSclzf2vxhvGgo3T/pC5xXW80Ln751ddzb1rOpZ
133                wzOXTp89l3iu1vic0WrTImOYugCd8uvrriy/aHf1w9nkizPPvDl/7rTe+cRTBmf3nVQx0T4DU0dMOAMA
134                p1CDeg==''')
135        with self.subTest('navier-stokes-pressure'):
136            self.assertAlmostEqual64(state1['p'], '''
137                eNoNiT0OQEAYBdd5xD1cSOIMiIJGVJuIiliS9VNYKolkNG6j9BUvmZlnmJnoSAnx6RmFNSUVjgErRVOQ
138                iFtWKTUZMQ2n/BuGnJaFi52bl48D73Fis+wjCpT6AWgsRHE=''')
139
140    def test_mixed(self):
141        state0, state1 = main(nelems=3, etype='mixed', reynolds=100, degree=2)
142        with self.subTest('stokes-velocity'):
143            self.assertAlmostEqual64(state0['u'], '''
144                eNpjYEAFEy++uHzs/EYjEFv73Hm9T2eVDGFyMWdfGJ/TVTGxMvQyqjxbAlYTZM7AsNWIxwhEG5hs0wTR
145                1Wd5DDTOHrx05OzmczC9Cud8riSccbrqYJR2dfMZ07M+hkznLpuongepB+Eyk/TzIHVbL4QbrLqw9Qyy