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