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
87        ax.add_collection(matplotlib.collections.LineCollection(x[bezier.hull], colors='w', linewidths=.5, alpha=.2))
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
120                VzP87Ce9b2f0DHaduW+w7EyEIYPeH4MW4z6Dh6apSOqKr4Wf5bv47cyl87vOKJ5fdmbFOQY9lvMtxkXn
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
146                m+aee31awWALXEzP+NIZkB6Y/QAD2Dbr''')
147        with self.subTest('stokes-pressure'):
148            self.assertAlmostEqual64(state0['p'], '''
149                eNoBIADf/1zNa81jzWHNs8w3zV3MUs2HzZrNzc26zanNd82qzgAAR/MTmQ==''')
150        with self.subTest('navier-stokes-velocity'):
151            self.assertAlmostEqual64(state1['u'], '''
152                eNpjYEAFvRcfXLa8eNsQxP5xLkBv17ktBjC5iHOfjbr1XxotNbAy0j4nbQQSCzJnYLA0lNIH0XLGi3VB
153                dNK5HoNjZ2frXj+77BxM775zfbq6Z+cYxBpd1vc543jursGSC0omz86D1IPwBJOzYDszL5oaLbgQcQbZ
154                TTZnec9svix4FsZXNbl9GqQHZj8AAcY1/g==''')
155        with self.subTest('navier-stokes-pressure'):
156            self.assertAlmostEqual64(state1['p'], '''
157                eNoBIADf/wXMDswMzAfMzMvry73LAcwIzCDM/ssJzDTM9MvRzAAAHqQRsw==''')