burgers.py

In this script we solve the Burgers equation on a 1D or 2D periodic domain, starting from a centered Gaussian and convecting in the positive direction of the first coordinate.

 7from nutils import mesh, function, solver, export, cli, testing
 8from nutils.expression_v2 import Namespace
 9import numpy
10import treelog
11
12
13def main(nelems: int, ndims: int, btype: str, degree: int, timescale: float, newtontol: float, endtime: float):
14    '''
15    Burgers equation on a 1D or 2D periodic domain.
16
17    .. arguments::
18
19       nelems [20]
20         Number of elements along a single dimension.
21       ndims [1]
22         Number of spatial dimensions.
23       btype [discont]
24         Type of basis function (discont/legendre).
25       degree [1]
26         Polynomial degree for discontinuous basis functions.
27       timescale [.5]
28         Fraction of timestep and element size: timestep=timescale/nelems.
29       newtontol [1e-5]
30         Newton tolerance.
31       endtime [inf]
32         Stopping time.
33    '''
34
35    domain, geom = mesh.rectilinear([numpy.linspace(-.5, .5, nelems+1)]*ndims, periodic=range(ndims))
36
37    ns = Namespace()
38    ns.x = geom
39    ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
40    ns.basis = domain.basis(btype, degree=degree)
41    ns.u = function.dotarg('lhs', ns.basis)
42    ns.f = '.5 u^2'
43    ns.C = 1
44    ns.u0 = 'exp(-25 x_i x_i)'
45
46    res = domain.integral('-∇_0(basis_n) f dV' @ ns, degree=5)
47    res += domain.interfaces.integral('-[basis_n] n_0 ({f} - .5 C [u] n_0) dS' @ ns, degree=degree*2)
48    inertia = domain.integral('basis_n u dV' @ ns, degree=5)
49
50    sqr = domain.integral('(u - u0)^2 dV' @ ns, degree=5)
51    lhs0 = solver.optimize('lhs', sqr)
52
53    timestep = timescale/nelems
54    bezier = domain.sample('bezier', 7)
55    with treelog.iter.plain('timestep', solver.impliciteuler('lhs', res, inertia, timestep=timestep, arguments=dict(lhs=lhs0), newtontol=newtontol)) as steps:
56        for itime, lhs in enumerate(steps):
57            x, u = bezier.eval(['x_i', 'u'] @ ns, lhs=lhs)
58            export.triplot('solution.png', x, u, tri=bezier.tri, hull=bezier.hull, clim=(0, 1))
59            if itime * timestep >= endtime:
60                break
61
62    return lhs

If the script is executed (as opposed to imported), nutils.cli.run() calls the main function with arguments provided from the command line. For example, to simulate until 0.5 seconds run python3 burgers.py endtime=0.5 (view log).

70if __name__ == '__main__':
71    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.

 80class test(testing.TestCase):
 81
 82    def test_1d_p0(self):
 83        lhs = main(ndims=1, nelems=10, timescale=.1, btype='discont', degree=0, endtime=.01, newtontol=1e-5)
 84        self.assertAlmostEqual64(lhs, '''
 85            eNrz1ttqGGOiZSZlrmbuZdZgcsEwUg8AOqwFug==''')
 86
 87    def test_1d_p1(self):
 88        lhs = main(ndims=1, nelems=10, timescale=.1, btype='discont', degree=1, endtime=.01, newtontol=1e-5)
 89        self.assertAlmostEqual64(lhs, '''
 90            eNrbocann6u3yqjTyMLUwfSw2TWzKPNM8+9mH8wyTMNNZxptMirW49ffpwYAI6cOVA==''')
 91
 92    def test_1d_p2(self):
 93        lhs = main(ndims=1, nelems=10, timescale=.1, btype='discont', degree=2, endtime=.01, newtontol=1e-5)
 94        self.assertAlmostEqual64(lhs, '''
 95            eNrr0c7SrtWfrD/d4JHRE6Ofxj6mnqaKZofNDpjZmQeYB5pHmL8we23mb5ZvWmjKY/LV6KPRFIMZ+o36
 96            8dp92gCxZxZG''')
 97
 98    def test_1d_p1_legendre(self):
 99        lhs = main(ndims=1, nelems=10, timescale=.1, btype='legendre', degree=1, endtime=.01, newtontol=1e-5)
100        self.assertAlmostEqual64(lhs, '''
101            eNrbpbtGt9VQyNDfxMdYzczERNZczdjYnOdsoNmc01kmE870Gj49t0c36BIAAhsO1g==''')
102
103    def test_1d_p2_legendre(self):
104        lhs = main(ndims=1, nelems=10, timescale=.1, btype='legendre', degree=2, endtime=.01, newtontol=1e-5)
105        self.assertAlmostEqual64(lhs, '''
106            eNoBPADD/8ot2y2/K4UxITFFLk00RTNNLyY2KzTTKx43QjOOzzM3Ss0pz1A2qsvhKGk0jsyXL48xzc5j
107            LswtIdLIK5SlF78=''')
108
109    def test_2d_p1(self):
110        lhs = main(ndims=2, nelems=4, timescale=.1, btype='discont', degree=1, endtime=.01, newtontol=1e-5)
111        import os
112        if os.environ.get('NUTILS_TENSORIAL'):
113            lhs = lhs.reshape(4, 2, 4, 2).transpose(0, 2, 1, 3).ravel()
114        self.assertAlmostEqual64(lhs, '''
115            eNoNyKENhEAQRuGEQsCv2SEzyQZHDbRACdsDJNsBjqBxSBxBHIgJ9xsqQJ1Drro1L1/eYBZceGz8njrR
116            yacm8UQLBvPYCw1airpyUVYSJLhKijK4IC01WDnqqxvX8OTl427aU73sctPGr3qqceBnRzOjo0xy9JpJ
117            R73m6R6YMZo/Q+FCLQ==''')