cylinderflow.py¶
In this script we solve the Navier-Stokes equations around a cylinder, using the same Raviart-Thomas discretization as in drivencavity-compatible.py but in curvilinear coordinates. The mesh is constructed such that all elements are shape similar, growing exponentially with radius such that the artificial exterior boundary is placed at large (configurable) distance.
10from nutils import mesh, function, solver, util, export, cli, testing
11from nutils.expression_v2 import Namespace
12import numpy
13import treelog
14
15
16def main(nelems: int, degree: int, reynolds: float, rotation: float, timestep: float, maxradius: float, seed: int, endtime: float):
17 '''
18 Flow around a cylinder.
19
20 .. arguments::
21
22 nelems [24]
23 Element size expressed in number of elements along the cylinder wall.
24 All elements have similar shape with approximately unit aspect ratio,
25 with elements away from the cylinder wall growing exponentially.
26 degree [3]
27 Polynomial degree for velocity space; the pressure space is one degree
28 less.
29 reynolds [1000]
30 Reynolds number, taking the cylinder radius as characteristic length.
31 rotation [0]
32 Cylinder rotation speed.
33 timestep [.04]
34 Time step
35 maxradius [25]
36 Target exterior radius; the actual domain size is subject to integer
37 multiples of the configured element size.
38 seed [0]
39 Random seed for small velocity noise in the intial condition.
40 endtime [inf]
41 Stopping time.
42 '''
43
44 elemangle = 2 * numpy.pi / nelems
45 melems = int(numpy.log(2*maxradius) / elemangle + .5)
46 treelog.info('creating {}x{} mesh, outer radius {:.2f}'.format(melems, nelems, .5*numpy.exp(elemangle*melems)))
47 domain, geom = mesh.rectilinear([melems, nelems], periodic=(1,))
48 domain = domain.withboundary(inner='left', outer='right')
49
50 ns = Namespace()
51 ns.δ = function.eye(domain.ndims)
52 ns.Σ = function.ones([domain.ndims])
53 ns.uinf = function.Array.cast([1, 0])
54 ns.r = .5 * function.exp(elemangle * geom[0])
55 ns.Re = reynolds
56 ns.phi = geom[1] * elemangle # add small angle to break element symmetry
57 ns.x_i = 'r (cos(phi) δ_i0 + sin(phi) δ_i1)'
58 ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
59 J = ns.x.grad(geom)
60 detJ = function.determinant(J)
61 ns.ubasis = function.matmat(function.vectorize([
62 domain.basis('spline', degree=(degree, degree-1), removedofs=((0,), None)),
63 domain.basis('spline', degree=(degree-1, degree))]), J.T) / detJ
64 ns.pbasis = domain.basis('spline', degree=degree-1) / detJ
65 ns.u = function.dotarg('u', ns.ubasis)
66 ns.p = function.dotarg('p', ns.pbasis)
67 ns.sigma_ij = '(∇_j(u_i) + ∇_i(u_j)) / Re - p δ_ij'
68 ns.N = 10 * degree / elemangle # Nitsche constant based on element size = elemangle/2
69 ns.nitsche_ni = '(N ubasis_ni - (∇_j(ubasis_ni) + ∇_i(ubasis_nj)) n_j) / Re'
70 ns.rotation = rotation
71 ns.uwall_i = '0.5 rotation (-sin(phi) δ_i0 + cos(phi) δ_i1)'
72
73 inflow = domain.boundary['outer'].select(-ns.uinf.dotnorm(ns.x), ischeme='gauss1') # upstream half of the exterior boundary
74 sqr = inflow.integral('Σ_i (u_i - uinf_i)^2' @ ns, degree=degree*2)
75 ucons = solver.optimize('u', sqr, droptol=1e-15) # constrain inflow semicircle to uinf
76 cons = dict(u=ucons)
77
78 numpy.random.seed(seed)
79 sqr = domain.integral('Σ_i (u_i - uinf_i)^2' @ ns, degree=degree*2)
80 udofs0 = solver.optimize('u', sqr) * numpy.random.normal(1, .1, len(ns.ubasis)) # set initial condition to u=uinf with small random noise
81 state0 = dict(u=udofs0)
82
83 ures = domain.integral('(ubasis_ni ∇_j(u_i) u_j + ∇_j(ubasis_ni) sigma_ij) dV' @ ns, degree=9)
84 ures += domain.boundary['inner'].integral('(nitsche_ni (u_i - uwall_i) - ubasis_ni sigma_ij n_j) dS' @ ns, degree=9)
85 pres = domain.integral('pbasis_n ∇_k(u_k) dV' @ ns, degree=9)
86 uinertia = domain.integral('ubasis_ni u_i dV' @ ns, degree=9)
87
88 bbox = numpy.array([[-2, 46/9], [-2, 2]]) # bounding box for figure based on 16x9 aspect ratio
89 bezier0 = domain.sample('bezier', 5)
90 bezier = bezier0.subset((bezier0.eval((ns.x-bbox[:, 0]) * (bbox[:, 1]-ns.x)) > 0).all(axis=1))
91 interpolate = util.tri_interpolator(bezier.tri, bezier.eval(ns.x), mergetol=1e-5) # interpolator for quivers
92 spacing = .05 # initial quiver spacing
93 xgrd = util.regularize(bbox, spacing)
94
95 with treelog.iter.plain('timestep', solver.impliciteuler(('u', 'p'), residual=(ures, pres), inertia=(uinertia, None), arguments=state0, timestep=timestep, constrain=cons, newtontol=1e-10)) as steps:
96 for istep, state in enumerate(steps):
97
98 t = istep * timestep
99 x, u, normu, p = bezier.eval(['x_i', 'u_i', 'sqrt(u_i u_i)', 'p'] @ ns, **state)
100 ugrd = interpolate[xgrd](u)
101
102 with export.mplfigure('flow.png', figsize=(12.8, 7.2)) as fig:
103 ax = fig.add_axes([0, 0, 1, 1], yticks=[], xticks=[], frame_on=False, xlim=bbox[0], ylim=bbox[1])
104 im = ax.tripcolor(x[:, 0], x[:, 1], bezier.tri, p, shading='gouraud', cmap='jet')
105 import matplotlib.collections
106 ax.add_collection(matplotlib.collections.LineCollection(x[bezier.hull], colors='k', linewidths=.1, alpha=.5))
107 ax.quiver(xgrd[:, 0], xgrd[:, 1], ugrd[:, 0], ugrd[:, 1], angles='xy', width=1e-3, headwidth=3e3, headlength=5e3, headaxislength=2e3, zorder=9, alpha=.5)
108 ax.plot(0, 0, 'k', marker=(3, 2, t*rotation*180/numpy.pi-90), markersize=20)
109 cax = fig.add_axes([0.9, 0.1, 0.01, 0.8])
110 cax.tick_params(labelsize='large')
111 fig.colorbar(im, cax=cax)
112
113 if t >= endtime:
114 break
115
116 xgrd = util.regularize(bbox, spacing, xgrd + ugrd * timestep)
117
118 return state0, state
If the script is executed (as opposed to imported), nutils.cli.run()
calls the main function with arguments provided from the command line.
124if __name__ == '__main__':
125 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.
134class test(testing.TestCase):
135
136 def test_rot0(self):
137 state0, state = main(nelems=6, degree=3, reynolds=100, rotation=0, timestep=.1, maxradius=25, seed=0, endtime=.05)
138 with self.subTest('initial condition'):
139 self.assertAlmostEqual64(state0['u'], '''
140 eNoNzLEJwkAYQOHXuYZgJyZHLjmMCIJY2WcEwQFcwcYFbBxCEaytzSVnuM7CQkEDqdVCG//uNe/rqEcI
141 p+pSwTzQAez90cM06kZQuLWDrV5oGJf9EupEJ7CyqYXUTAyM7C2HmZyNf8p57S2lB95It8KNgmH1PcNB
142 /UTEvUVpojqGdpEV8IlfIt7znSiZ+QPSaDIR''', atol=2e-13)
143 with self.subTest('velocity'):
144 self.assertAlmostEqual64(state['u'], '''
145 eNoBkABv/+o0szWg04bKlsogMVI4JjcmMXXI+cfizb05/Dk4MBHGEcaPNDo8ljuTNibE4sNpznI9VD02
146 M5zCnsJazE0+Hj76NsPByMH/yl43nDNlyGnIsy+YNz44MspbyD/IWcwHOMM4AzXAxsPGGjAJOUM7GcgU
147 xePEqckvO+g8+DcOwyfD4zfjPFY+sMfJwavBhDNPPpLTRUE=''')
148 with self.subTest('pressure'):
149 self.assertAlmostEqual64(state['p'], '''
150 eNoBSAC3/4zF18aozR866DpHNSk8JDonOdw4k8VzNaHBk8PFOyI+Gj9vPPRA/T/LQDtBIECaP0i5yLsA
151 wL9FwkabQsJJTbc2ubJHw7ZvRq/qITA=''')
152
153 def test_rot1(self):
154 state0, state = main(nelems=6, degree=3, reynolds=100, rotation=1, timestep=.1, maxradius=25, seed=0, endtime=.05)
155 with self.subTest('initial condition'):
156 self.assertAlmostEqual64(state0['u'], '''
157 eNoNzLEJwkAYQOHXuYZgJyZHLjmMCIJY2WcEwQFcwcYFbBxCEaytzSVnuM7CQkEDqdVCG//uNe/rqEcI
158 p+pSwTzQAez90cM06kZQuLWDrV5oGJf9EupEJ7CyqYXUTAyM7C2HmZyNf8p57S2lB95It8KNgmH1PcNB
159 /UTEvUVpojqGdpEV8IlfIt7znSiZ+QPSaDIR''', atol=2e-13)
160 with self.subTest('velocity'):
161 self.assertAlmostEqual64(state['u'], '''
162 eNoBkABv/+M0tzVcKYfKlMr1MFE4JzdBMXbI+cfMzb05/Dk/MBHGEcaPNDo8ljuUNibE4sNjznI9VD02
163 M5zCnsJazE0+Hj76NsPByMH/ynk3WTR/yH7I5TGsNzk4HcpUyDnILMwBOMQ4BzXBxsTGpDAKOUM7GMgU
164 xePEp8kvO+g8+DcOwyfD5DfkPFY+sMfJwavBhDNPPpo5RbI=''')
165 with self.subTest('pressure'):
166 self.assertAlmostEqual64(state['p'], '''
167 eNoBSAC3/4rF3Ma4ziE65zoUNSg8JToqOd04k8VbNaDBlMPJOyI+Gj9sPPRA/T/MQDtBH0CYP0e5yLsD
168 wL9FwkaaQsJJTbc3ubJHw7ZvRqV7IP8=''')