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
 11import numpy, treelog
 12
 13def main(nelems:int, degree:int, reynolds:float, rotation:float, timestep:float, maxradius:float, seed:int, endtime:float):
 14  '''
 15  Flow around a cylinder.
 16
 17  .. arguments::
 18
 19     nelems [24]
 20       Element size expressed in number of elements along the cylinder wall.
 21       All elements have similar shape with approximately unit aspect ratio,
 22       with elements away from the cylinder wall growing exponentially.
 23     degree [3]
 24       Polynomial degree for velocity space; the pressure space is one degree
 25       less.
 26     reynolds [1000]
 27       Reynolds number, taking the cylinder radius as characteristic length.
 28     rotation [0]
 29       Cylinder rotation speed.
 30     timestep [.04]
 31       Time step
 32     maxradius [25]
 33       Target exterior radius; the actual domain size is subject to integer
 34       multiples of the configured element size.
 35     seed [0]
 36       Random seed for small velocity noise in the intial condition.
 37     endtime [inf]
 38       Stopping time.
 39  '''
 40
 41  elemangle = 2 * numpy.pi / nelems
 42  melems = int(numpy.log(2*maxradius) / elemangle + .5)
 43  treelog.info('creating {}x{} mesh, outer radius {:.2f}'.format(melems, nelems, .5*numpy.exp(elemangle*melems)))
 44  domain, geom = mesh.rectilinear([melems, nelems], periodic=(1,))
 45  domain = domain.withboundary(inner='left', outer='right')
 46
 47  ns = function.Namespace()
 48  ns.uinf = 1, 0
 49  ns.r = .5 * function.exp(elemangle * geom[0])
 50  ns.Re = reynolds
 51  ns.phi = geom[1] * elemangle # add small angle to break element symmetry
 52  ns.x_i = 'r <cos(phi), sin(phi)>_i'
 53  ns.J = ns.x.grad(geom)
 54  ns.unbasis, ns.utbasis, ns.pbasis = function.chain([ # compatible spaces
 55    domain.basis('spline', degree=(degree,degree-1), removedofs=((0,),None)),
 56    domain.basis('spline', degree=(degree-1,degree)),
 57    domain.basis('spline', degree=degree-1),
 58  ]) / function.determinant(ns.J)
 59  ns.ubasis_ni = 'unbasis_n J_i0 + utbasis_n J_i1' # piola transformation
 60  ns.u_i = 'ubasis_ni ?lhs_n'
 61  ns.p = 'pbasis_n ?lhs_n'
 62  ns.sigma_ij = '(u_i,j + u_j,i) / Re - p δ_ij'
 63  ns.N = 10 * degree / elemangle # Nitsche constant based on element size = elemangle/2
 64  ns.nitsche_ni = '(N ubasis_ni - (ubasis_ni,j + ubasis_nj,i) n_j) / Re'
 65  ns.rotation = rotation
 66  ns.uwall_i = '0.5 rotation <-sin(phi), cos(phi)>_i'
 67
 68  inflow = domain.boundary['outer'].select(-ns.uinf.dotnorm(ns.x), ischeme='gauss1') # upstream half of the exterior boundary
 69  sqr = inflow.integral('(u_i - uinf_i) (u_i - uinf_i)' @ ns, degree=degree*2)
 70  cons = solver.optimize('lhs', sqr, droptol=1e-15) # constrain inflow semicircle to uinf
 71
 72  sqr = domain.integral('(u_i - uinf_i) (u_i - uinf_i) + p^2' @ ns, degree=degree*2)
 73  lhs0 = solver.optimize('lhs', sqr) # set initial condition to u=uinf, p=0
 74
 75  numpy.random.seed(seed)
 76  lhs0 *= numpy.random.normal(1, .1, lhs0.shape) # add small velocity noise
 77
 78  res = domain.integral('(ubasis_ni u_i,j u_j + ubasis_ni,j sigma_ij + pbasis_n u_k,k) d:x' @ ns, degree=9)
 79  res += domain.boundary['inner'].integral('(nitsche_ni (u_i - uwall_i) - ubasis_ni sigma_ij n_j) d:x' @ ns, degree=9)
 80  inertia = domain.integral('ubasis_ni u_i d:x' @ ns, degree=9)
 81
 82  bbox = numpy.array([[-2,46/9],[-2,2]]) # bounding box for figure based on 16x9 aspect ratio
 83  bezier0 = domain.sample('bezier', 5)
 84  bezier = bezier0.subset((bezier0.eval((ns.x-bbox[:,0]) * (bbox[:,1]-ns.x)) > 0).all(axis=1))
 85  interpolate = util.tri_interpolator(bezier.tri, bezier.eval(ns.x), mergetol=1e-5) # interpolator for quivers
 86  spacing = .05 # initial quiver spacing
 87  xgrd = util.regularize(bbox, spacing)
 88
 89  with treelog.iter.plain('timestep', solver.impliciteuler('lhs', residual=res, inertia=inertia, lhs0=lhs0, timestep=timestep, constrain=cons, newtontol=1e-10)) as steps:
 90    for istep, lhs in enumerate(steps):
 91
 92      t = istep * timestep
 93      x, u, normu, p = bezier.eval(['x_i', 'u_i', 'sqrt(u_k u_k)', 'p'] @ ns, lhs=lhs)
 94      ugrd = interpolate[xgrd](u)
 95
 96      with export.mplfigure('flow.png', figsize=(12.8,7.2)) as fig:
 97        ax = fig.add_axes([0,0,1,1], yticks=[], xticks=[], frame_on=False, xlim=bbox[0], ylim=bbox[1])
 98        im = ax.tripcolor(x[:,0], x[:,1], bezier.tri, p, shading='gouraud', cmap='jet')
 99        import matplotlib.collections
100        ax.add_collection(matplotlib.collections.LineCollection(x[bezier.hull], colors='k', linewidths=.1, alpha=.5))
101        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)
102        ax.plot(0, 0, 'k', marker=(3,2,t*rotation*180/numpy.pi-90), markersize=20)
103        cax = fig.add_axes([0.9, 0.1, 0.01, 0.8])
104        cax.tick_params(labelsize='large')
105        fig.colorbar(im, cax=cax)
106
107      if t >= endtime:
108        break
109
110      xgrd = util.regularize(bbox, spacing, xgrd + ugrd * timestep)
111
112  return lhs0, lhs

If the script is executed (as opposed to imported), nutils.cli.run() calls the main function with arguments provided from the command line.

117if __name__ == '__main__':
118  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.

126class test(testing.TestCase):
127
128  @testing.requires('matplotlib', 'scipy')
129  def test_rot0(self):
130    lhs0, lhs = main(nelems=6, degree=3, reynolds=100, rotation=0, timestep=.1, maxradius=25, seed=0, endtime=.05)
131    with self.subTest('initial condition'): self.assertAlmostEqual64(lhs0, '''
132      eNpT1n+qx8Bw8sLNCwwM6bpGugwMmy7tv8TA4GmoZcjAcObctHMMDOuNio0YGBzPmp9lYHhuYmTCwNB5
133      2uI0A4OFqbMpA4Pd6YenGBhSgDpfXXoG1HlXpwXItrxkCmSz683WZ2CwvvDrPAPDVv3fQBMZzn0FmvLK
134      8LkxA4PCmZAzDAzfjL8ATXx0agPQlBCgedQBAOgCMhE=''', atol=2e-13)
135    with self.subTest('left-hand side'): self.assertAlmostEqual64(lhs, '''
136      eNoB2AAn/4Y0pjUwMHTKhMrmMoI4Qzcpz4TI78egy545+Dm7MwPGEsa+NVY8pjtVNSzE18OoyXI9VD02
137      M5zCnsJazE0+Hj76NsPByMH/yhQ30DN6yFjIAjCrN5Y4FcooyE3I8ssCOGk4QjXXxrPGNzILOXo7AMj3
138      xOjEM8k3O8Y85DcZwyTDAzjaPFY+sMfJwavBhDNPPvbFX8cuOKI3/zpFOFI87TqmN9k8C8hkNFnCgcXV
139      Pds7VT/qPdZBbEF5QUZD7UEJQYi527ziROVETEeVRfZIfrfuRKZKr7s6SRCVaAA=''')
140
141  @testing.requires('matplotlib', 'scipy')
142  def test_rot1(self):
143    lhs0, lhs = main(nelems=6, degree=3, reynolds=100, rotation=1, timestep=.1, maxradius=25, seed=0, endtime=.05)
144    with self.subTest('initial condition'): self.assertAlmostEqual64(lhs0, '''
145      eNpT1n+qx8Bw8sLNCwwM6bpGugwMmy7tv8TA4GmoZcjAcObctHMMDOuNio0YGBzPmp9lYHhuYmTCwNB5
146      2uI0A4OFqbMpA4Pd6YenGBhSgDpfXXoG1HlXpwXItrxkCmSz683WZ2CwvvDrPAPDVv3fQBMZzn0FmvLK
147      8LkxA4PCmZAzDAzfjL8ATXx0agPQlBCgedQBAOgCMhE=''', atol=2e-13)
148    with self.subTest('left-hand side'): self.assertAlmostEqual64(lhs, '''
149      eNoB2AAn/380qTWFMHXKgsrUMoI4RDdJz4XI78eZy545+Dm8MwPGEsa+NVY8pjtWNSzE18OoyXI9VD02
150      M5zCnsJazE0+Hj76NsPByMH/yjM3ejSWyGzI/TG+N5I4A8oiyEjIzsv9N2o4RTXYxrTGajIMOXo7AMj3
151      xOjEMsk3O8Y85TcZwyTDAzjaPFY+sMfJwavBhDNPPvPFZMc4OKg3/jo7OFI87jqtN9k8Ccg6NFjChcXW
152      Pd07VT/oPdZBbEF5QUZD7EEIQYe527ziROVETEeURfZIfrfuRKZKrrs6SVFLajU=''')