# 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, '''