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
13import numpy, treelog
14
15def main(nelems:int, degree:int, reynolds:float):
16 '''
17 Driven cavity benchmark problem using compatible spaces.
18
19 .. arguments::
20
21 nelems [12]
22 Number of elements along edge.
23 degree [2]
24 Polynomial degree for velocity; the pressure space is one degree less.
25 reynolds [1000]
26 Reynolds number, taking the domain size as characteristic length.
27 '''
28
29 verts = numpy.linspace(0, 1, nelems+1)
30 domain, geom = mesh.rectilinear([verts, verts])
31
32 ns = function.Namespace()
33 ns.x = geom
34 ns.Re = reynolds
35 ns.uxbasis, ns.uybasis, ns.pbasis, ns.lbasis = function.chain([
36 domain.basis('spline', degree=(degree,degree-1), removedofs=((0,-1),None)),
37 domain.basis('spline', degree=(degree-1,degree), removedofs=(None,(0,-1))),
38 domain.basis('spline', degree=degree-1),
39 [1], # lagrange multiplier
40 ])
41 ns.ubasis_ni = '<uxbasis_n, uybasis_n>_i'
42 ns.u_i = 'ubasis_ni ?lhs_n'
43 ns.p = 'pbasis_n ?lhs_n'
44 ns.l = 'lbasis_n ?lhs_n'
45 ns.stress_ij = '(u_i,j + u_j,i) / Re - p δ_ij'
46 ns.uwall = domain.boundary.indicator('top'), 0
47 ns.N = 5 * degree * nelems # nitsche constant based on element size = 1/nelems
48 ns.nitsche_ni = '(N ubasis_ni - (ubasis_ni,j + ubasis_nj,i) n_j) / Re'
49
50 res = domain.integral('(ubasis_ni,j stress_ij + pbasis_n (u_k,k + l) + lbasis_n p) d:x' @ ns, degree=2*degree)
51 res += domain.boundary.integral('(nitsche_ni (u_i - uwall_i) - ubasis_ni stress_ij n_j) d:x' @ ns, degree=2*degree)
52 with treelog.context('stokes'):
53 lhs0 = solver.solve_linear('lhs', res)
54 postprocess(domain, ns, lhs=lhs0)
55
56 res += domain.integral('ubasis_ni u_i,j u_j d:x' @ ns, degree=3*degree)
57 with treelog.context('navierstokes'):
58 lhs1 = solver.newton('lhs', res, lhs0=lhs0).solve(tol=1e-10)
59 postprocess(domain, ns, lhs=lhs1)
60
61 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.
67def postprocess(domain, ns, every=.05, spacing=.01, **arguments):
68
69 div = domain.integral('(u_k,k)^2 d:x' @ ns, degree=1).eval(**arguments)**.5
70 treelog.info('velocity divergence: {:.2e}'.format(div)) # confirm that velocity is pointwise divergence-free
71
72 ns = ns.copy_() # copy namespace so that we don't modify the calling argument
73 ns.streambasis = domain.basis('std', degree=2)[1:] # remove first dof to obtain non-singular system
74 ns.stream = 'streambasis_n ?streamdofs_n' # stream function
75 sqr = domain.integral('((u_0 - stream_,1)^2 + (u_1 + stream_,0)^2) d:x' @ ns, degree=4)
76 arguments['streamdofs'] = solver.optimize('streamdofs', sqr, arguments=arguments) # compute streamlines
77
78 bezier = domain.sample('bezier', 9)
79 x, u, p, stream = bezier.eval(['x_i', 'sqrt(u_i u_i)', 'p', 'stream'] @ ns, **arguments)
80 with export.mplfigure('flow.png') as fig: # plot velocity as field, pressure as contours, streamlines as dashed
81 ax = fig.add_axes([.1,.1,.8,.8], yticks=[], aspect='equal')
82 import matplotlib.collections
83 ax.add_collection(matplotlib.collections.LineCollection(x[bezier.hull], colors='w', linewidths=.5, alpha=.2))
84 ax.tricontour(x[:,0], x[:,1], bezier.tri, stream, 16, colors='k', linestyles='dotted', linewidths=.5, zorder=9)
85 caxu = fig.add_axes([.1,.1,.03,.8], title='velocity')
86 imu = ax.tripcolor(x[:,0], x[:,1], bezier.tri, u, shading='gouraud', cmap='jet')
87 fig.colorbar(imu, cax=caxu)
88 caxu.yaxis.set_ticks_position('left')
89 caxp = fig.add_axes([.87,.1,.03,.8], title='pressure')
90 imp = ax.tricontour(x[:,0], x[:,1], bezier.tri, p, 16, cmap='gray', linestyles='solid')
91 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).
98if __name__ == '__main__':
99 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.
107class test(testing.TestCase):
108
109 @testing.requires('matplotlib')
110 def test_p1(self):
111 lhs0, lhs1 = main(nelems=3, reynolds=100, degree=2)
112 with self.subTest('stokes'): self.assertAlmostEqual64(lhs0, '''
113 eNrzu9Bt8OuUndkD/eTTSqezzP2g/E3698/ZmZlf2GjSaHJS3/90/Wm/C4qGh066XzLQ47846VSPpoWK
114 3vnD+iXXTty+ZGB7YafuhYsf9fJMGRgAkFIn4A==''')
115 with self.subTest('navier-stokes'): self.assertAlmostEqual64(lhs1, '''
116 eNoBUgCt/2XOWjJSy5k1jS+yyzvLODfgL1rO0MrINpsxHM2ZNSrPqDTANCPVQsxCzeAvcc04yaUmYysm
117 MbLLAi9YL6TN+y3eLcgvM87NzOUrTNY9MWA2AABnnyYn''')
118
119 @testing.requires('matplotlib')
120 def test_p2(self):
121 lhs0, lhs1 = main(nelems=3, reynolds=100, degree=3)
122 with self.subTest('stokes'): self.assertAlmostEqual64(lhs0, '''
123 eNo7aLjtjIjJxZN7zVgvZJ9jOv3lfK05gnUQLmt/Ttlk5qm9ZgKGQeeXmj0zZoCCD+fWGUSflDpz0PDu
124 6XRT55OL9dt11pwvNYw5+f7ClYv2Oq/O7DBigANBfR29g5fFjD3Oxl6ovBxi0H1uiRkDAwD+ITkl''')
125 with self.subTest('navier-stokes'): self.assertAlmostEqual64(lhs1, '''
126 eNoBhAB7/14yGcxyNPbJYTahLj/LSDE7yy43SM9WMsXJoDR+N3Iw8s1hM5zJODeizcE0X8phNrQwUDOO
127 NbMzJi+ty4s1oDFqzxIzysjWzXIwFM3tNMjIKjG8MeLNoTLzyQMuCi+IK3jOcjMzLuMvudNEzOrOEDAF
128 MD8sTTDpzNjYZDCg0RgwcTcAAJCyOzM=''')