drivencavity.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.
8 | import nutils, numpy
|
The main function defines the parameter space for the script. Configurable parameters are the mesh density (in number of elements along an edge), element type (square, triangle, or mixed), polynomial degree, and Reynolds number.
15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 | def main(nelems: 'number of elements' = 12,
etype: 'type of elements (square/triangle/mixed)' = 'square',
degree: 'polynomial degree for velocity' = 3,
reynolds: 'reynolds number' = 1000.):
domain, geom = nutils.mesh.unitsquare(nelems, etype)
ns = nutils.function.Namespace()
ns.Re = reynolds
ns.x = geom
ns.ubasis, ns.pbasis = nutils.function.chain([
domain.basis('std', degree=degree).vector(2),
domain.basis('std', degree=degree-1),
])
ns.u_i = 'ubasis_ni ?lhs_n'
ns.p = 'pbasis_n ?lhs_n'
ns.stress_ij = '(u_i,j + u_j,i) / Re - p δ_ij'
sqr = domain.boundary.integral('u_k u_k d:x' @ ns, degree=degree*2)
wallcons = nutils.solver.optimize('lhs', sqr, droptol=1e-15)
sqr = domain.boundary['top'].integral('(u_0 - 1)^2 d:x' @ ns, degree=degree*2)
lidcons = nutils.solver.optimize('lhs', sqr, droptol=1e-15)
cons = numpy.choose(numpy.isnan(lidcons), [lidcons, wallcons])
cons[-1] = 0 # pressure point constraint
res = domain.integral('(ubasis_ni,j stress_ij + pbasis_n u_k,k) d:x' @ ns, degree=degree*2)
with nutils.log.context('stokes'):
lhs0 = nutils.solver.solve_linear('lhs', res, constrain=cons)
postprocess(domain, ns, lhs=lhs0)
res += domain.integral('ubasis_ni u_i,j u_j d:x' @ ns, degree=degree*3)
with nutils.log.context('navierstokes'):
lhs1 = nutils.solver.newton('lhs', res, lhs0=lhs0, constrain=cons).solve(tol=1e-10)
postprocess(domain, ns, lhs=lhs1)
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.
58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 | def postprocess(domain, ns, every=.05, spacing=.01, **arguments):
ns = ns.copy_() # copy namespace so that we don't modify the calling argument
ns.streambasis = domain.basis('std', degree=2)[1:] # remove first dof to obtain non-singular system
ns.stream = 'streambasis_n ?streamdofs_n' # stream function
sqr = domain.integral('((u_0 - stream_,1)^2 + (u_1 + stream_,0)^2) d:x' @ ns, degree=4)
arguments['streamdofs'] = nutils.solver.optimize('streamdofs', sqr, arguments=arguments) # compute streamlines
bezier = domain.sample('bezier', 9)
x, u, p, stream = bezier.eval(['x_i', 'sqrt(u_k u_k)', 'p', 'stream'] @ ns, **arguments)
with nutils.export.mplfigure('flow.png') as fig: # plot velocity as field, pressure as contours, streamlines as dashed
ax = fig.add_axes([.1,.1,.8,.8], yticks=[], aspect='equal')
import matplotlib.collections
ax.add_collection(matplotlib.collections.LineCollection(x[bezier.hull], colors='w', linewidths=.5, alpha=.2))
ax.tricontour(x[:,0], x[:,1], bezier.tri, stream, 16, colors='k', linestyles='dotted', linewidths=.5, zorder=9)
caxu = fig.add_axes([.1,.1,.03,.8], title='velocity')
imu = ax.tripcolor(x[:,0], x[:,1], bezier.tri, u, shading='gouraud', cmap='jet')
fig.colorbar(imu, cax=caxu)
caxu.yaxis.set_ticks_position('left')
caxp = fig.add_axes([.87,.1,.03,.8], title='pressure')
imp = ax.tricontour(x[:,0], x[:,1], bezier.tri, p, 16, cmap='gray', linestyles='solid')
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.py
(view log).
85 86 | if __name__ == '__main__':
nutils.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.
94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 | class test(nutils.testing.TestCase):
@nutils.testing.requires('matplotlib')
def test_square(self):
lhs0, lhs1 = main(nelems=3, etype='square', reynolds=100, degree=3)
with self.subTest('stokes'): self.assertAlmostEqual64(lhs0, '''
eNp1zj1IQlEUB/BrCJKEQxLRFNFQxvN1vTcpoqWhzZaGElr7WKOGirApiIaipcEKoiXCpaKEiCKnhjzn
XX1PejaEJGGFRCCiCH153YrXOXCG3+Fw/oT8rZFeQpaVqDGVmjHNxEKSJmxM2rOIal1aDlsxKyK+gF/a
sZbHEA5gDmL6FduuWRnHsAQXcABEXeGP/5rVrdUPqyxWma1q2ih3u1g7/+JnPf3+BiYtr5ToBGvm33yN
d/C3pLTrTi9d9Y2yCkuxU2Z6pa17CqpKMzTo+6AbdLJmc3eupC7axKFmF7NiR5c2aBpiUYugAxUcRk/N
mgyn2MVXsME83INblRZW6hMFfIA6CMRvbotonTgL7/ACWQjBfjwcT8MT6HAJSxCEI8hAvroxIQZ7cA7F
X+3ET3CgG1Ucxz5sRDu2IMctTONQNVkFbNW5iScGIT8HbdXq''')
with self.subTest('navier-stokes'): self.assertAlmostEqual64(lhs1, '''
eNptzktoU0EUBuC7KeLGguKioS4MBdPekNyZSWIwEihowVVBxJW0pYuiFgpiXSh0F0ltELvoC2zAVuor
RuiTJlRLC6Hof2cml0wwCxVqCl1XFOqi4p27LPlXP985HI5hHM/1i4aRMzvVL7VqOs4j5VMhS9un8k2Z
kEnZLL+271v3mLYb8oG4KuKiR0yGtkk6om1MODzLH/Ma/xZK0b+eXROveJzX7Vs8ZcXYUFTbkYiJp7yF
b9i3VTO765m/fFL+5IM8ZBfFHJvybCD4WvVWi86BZPIsj3j3Gv3cKKXKUDhJovQ7TbBhdsrSdjl4xcqS
btrEZukM7VDa3ge2wnHSRAt0lmboSFjbCfNMuGItkH7aSxdpi9Q2c+Gf80JFgpdIHxkgdaJtt3aufFq2
iRXxUPqchLfnV63yLT/Pd2CKLXqfadsL9DmGmLeruPPl42diN/44jyV8wBuMogvteIe827MYxwTWkMOi
K1k8QxrTbl9xZQpPMIzn2EDR3cgjg5dYxzYKKIHjDzbx252sY9mdHuKHaRj/AYh1yFc=''')
@nutils.testing.requires('matplotlib')
def test_mixed(self):
lhs0, lhs1 = main(nelems=3, etype='mixed', reynolds=100, degree=2)
with self.subTest('stokes'): self.assertAlmostEqual64(lhs0, '''
eNpjYICAiRePnWdg0D736SyIF3P2nK6VYSWQHWS+1SjI3MAkyLz6rMbZI2BZhXMJZxyMNp/xMbwMFA8y
LzNhYNh6YdUFiElzzykYgGg94yBzkH6oBQwvLm80YmA4r6dkCOYZq5h4GZUYgdg8QHKbJpA2OHhp8zmQ
iM8Vp6tpV03PMp1TPQ/ipwPJcIOtZyAmvT69Bcy6BOXHnM0+m3w28ezmM+ZnY88EnW0/O+vs2bO7zq48
W352FdA8ABC3SoM=''')
with self.subTest('navier-stokes'): self.assertAlmostEqual64(lhs1, '''
eNpjYICA1RezLjIwPD639hyIl31umX6vgQGQHWTuaRhkLmYcZB54bvvZq2dBsofPqZ4tMoo4o22oaxJk
HmReasLAsOrihAsQkxzOJl0B0TJAOZB+qAUMtZefGzIwxOjtNgDxfho9MbI1UjcCsV/pMTA802VgqDNY
qrsEbL+I7nGD0/o655ouMIFN3QLUqWSUcQZiEvMZbrA7npyG8IXPyJ2RPiN65ubpn6dPn+Y9I3XG4Awf
UMzlDPuZ60A9AH73RT0=''')
|