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
91 ax.add_collection(matplotlib.collections.LineCollection(x[bezier.hull], colors='w', linewidths=.5, alpha=.2))
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'], '''
126 eNoBIADf/0fSMC4P0ZLKjCk4JC7Pwy901sjb0jA90Lkt0NHxLm41+KgP+Q==''')
127 with self.subTest('stokes-multiplier'):
128 self.assertAlmostEqual(state0['lm'], 0)
129 with self.subTest('navier-stokes-velocity'):
130 self.assertAlmostEqual64(state1['u'], '''
131 eNoBMADP/2XOWjJSy5k1jS+yyzvLODfgL1rO0MrINpsxHM2ZNSrPqDTANCPVQsxCzeAvcc04yT3AF9c=''')
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)