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