Merge pull request #682 from evalf/fields

Fields

examples/adaptivity.py

@@ -42,44 +42,45 @@ def main(etype: str, btype: str, degree: int, nrefine: int):
     domain = domain.trim(exact-1e-15, maxrefine=0)
     linreg = util.linear_regressor()
 
-    with treelog.iter.fraction('level', range(nrefine+1)) as lrange:
-        for irefine in lrange:
-
-            if irefine:
-                refdom = domain.refined
-                ns.refbasis = refdom.basis(btype, degree=degree)
-                indicator = refdom.integral('∇_k(refbasis_n) ∇_k(u) dV' @ ns, degree=degree*2).eval(lhs=lhs)
-                indicator -= refdom.boundary.integral('refbasis_n ∇_k(u) n_k dS' @ ns, degree=degree*2).eval(lhs=lhs)
-                supp = ns.refbasis.get_support(indicator**2 > numpy.mean(indicator**2))
-                domain = domain.refined_by(refdom.transforms[supp])
-
-            ns = Namespace()
-            ns.x = geom
-            ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
-            ns.basis = domain.basis(btype, degree=degree)
-            ns.u = function.dotarg('lhs', ns.basis)
-            ns.du = ns.u - exact
-
-            sqr = domain.boundary['trimmed'].integral('u^2 dS' @ ns, degree=degree*2)
-            cons = solver.optimize('lhs', sqr, droptol=1e-15)
-
-            sqr = domain.boundary.integral('du^2 dS' @ ns, degree=7)
-            cons = solver.optimize('lhs', sqr, droptol=1e-15, constrain=cons)
-
-            res = domain.integral('∇_k(basis_n) ∇_k(u) dV' @ ns, degree=degree*2)
-            lhs = solver.solve_linear('lhs', res, constrain=cons)
-
-            ndofs = len(ns.basis)
-            error = function.sqrt(domain.integral(['du du dV', '∇_k(du) ∇_k(du) dV'] @ ns, degree=7)).eval(lhs=lhs)
-            rate, offset = linreg.add(numpy.log(len(ns.basis)), numpy.log(error))
-            treelog.user('ndofs: {ndofs}, L2 error: {error[0]:.2e} ({rate[0]:.2f}), H1 error: {error[1]:.2e} ({rate[1]:.2f})'.format(ndofs=len(ns.basis), error=error, rate=rate))
-
-            bezier = domain.sample('bezier', 9)
-            x, u, du = bezier.eval(['x_i', 'u', 'du'] @ ns, lhs=lhs)
-            export.triplot('sol.png', x, u, tri=bezier.tri, hull=bezier.hull)
-            export.triplot('err.png', x, du, tri=bezier.tri, hull=bezier.hull)
-
-    return ndofs, error, lhs
+    for irefine in treelog.iter.fraction('level', range(nrefine+1)):
+
+        if irefine:
+            refdom = domain.refined
+            refbasis = refdom.basis(btype, degree=degree)
+            ns.add_field('vref', refbasis)
+            res = refdom.integral('∇_k(vref) ∇_k(u) dV' @ ns, degree=degree*2)
+            res -= refdom.boundary.integral('vref ∇_k(u) n_k dS' @ ns, degree=degree*2)
+            indicator = res.derivative('vref').eval(**args)
+            supp = refbasis.get_support(indicator**2 > numpy.mean(indicator**2))
+            domain = domain.refined_by(refdom.transforms[supp])
+
+        ns = Namespace()
+        ns.x = geom
+        ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
+        ns.add_field(('u', 'v'), domain.basis(btype, degree=degree))
+        ns.uexact = exact
+        ns.du = 'u - uexact'
+
+        sqr = domain.boundary['trimmed'].integral('u^2 dS' @ ns, degree=degree*2)
+        cons = solver.optimize('u,', sqr, droptol=1e-15)
+
+        sqr = domain.boundary.integral('du^2 dS' @ ns, degree=7)
+        cons = solver.optimize('u,', sqr, droptol=1e-15, constrain=cons)
+
+        res = domain.integral('∇_k(v) ∇_k(u) dV' @ ns, degree=degree*2)
+        args = solver.solve_linear('u:v', res, constrain=cons)
+
+        ndofs = len(args['u'])
+        error = function.sqrt(domain.integral(['du du dV', '∇_k(du) ∇_k(du) dV'] @ ns, degree=7)).eval(**args)
+        rate, offset = linreg.add(numpy.log(ndofs), numpy.log(error))
+        treelog.user(f'ndofs: {ndofs}, L2 error: {error[0]:.2e} ({rate[0]:.2f}), H1 error: {error[1]:.2e} ({rate[1]:.2f})')
+
+        bezier = domain.sample('bezier', 9)
+        xsmp, usmp, dusmp = bezier.eval(['x_i', 'u', 'du'] @ ns, **args)
+        export.triplot('sol.png', xsmp, usmp, tri=bezier.tri, hull=bezier.hull)
+        export.triplot('err.png', xsmp, dusmp, tri=bezier.tri, hull=bezier.hull)
+
+    return error, args['u']
 
 # If the script is executed (as opposed to imported), :func:`nutils.cli.run`
 # calls the main function with arguments provided from the command line. For
@@ -101,44 +102,44 @@ def main(etype: str, btype: str, degree: int, nrefine: int):
 class test(testing.TestCase):
 
     def test_square_quadratic(self):
-        ndofs, error, lhs = main(nrefine=2, btype='h-std', etype='square', degree=2)
+        error, u = main(nrefine=2, btype='h-std', etype='square', degree=2)
         with self.subTest('degrees of freedom'):
-            self.assertEqual(ndofs, 149)
+            self.assertEqual(len(u), 149)
         with self.subTest('L2-error'):
             self.assertAlmostEqual(error[0], 0.00065, places=5)
         with self.subTest('H1-error'):
             self.assertAlmostEqual(error[1], 0.03461, places=5)
         with self.subTest('left-hand side'):
-            self.assertAlmostEqual64(lhs, '''
+            self.assertAlmostEqual64(u, '''
                 eNo1j6FrQmEUxT8RBi4KllVfMsl3z/nK4zEmLC6bhsKCw2gSw5IPFsymGbZiWnr+By8Ii7Yhsk3BMtC4
                 Z9sJ223ncs85vzvmM9+Yhix8hDIjtnkdHqQSdDDDj1Qajr5qPXN/07MZ2vI4V7UOIvmdO/oEZY45xYDn
                 oR7ikLHAHVpcs2A1TLhChDO+MOeWt5xjYzm6fOQrGxxiZPeoMGaf37hCyU72hB0u6PglPcQcKxRI/KUd
                 7AYLvMPpsqGkCTPumzWf+qV92kKevjK36ozDP/FSnh1iteWiqWuf+oMaKuyKaC1i52rKPokiF2WLA/20
                 bya+ZCPbWKRPpvgFaedebw==''')
 
     def test_triangle_quadratic(self):
-        ndofs, error, lhs = main(nrefine=2, btype='h-std', etype='triangle', degree=2)
+        error, u = main(nrefine=2, btype='h-std', etype='triangle', degree=2)
         with self.subTest('degrees of freedom'):
-            self.assertEqual(ndofs, 98)
+            self.assertEqual(len(u), 98)
         with self.subTest('L2-error'):
             self.assertAlmostEqual(error[0], 0.00138, places=5)
         with self.subTest('H1-error'):
             self.assertAlmostEqual(error[1], 0.05324, places=5)
         with self.subTest('left-hand side'):
-            self.assertAlmostEqual64(lhs, '''
+            self.assertAlmostEqual64(u, '''
                 eNprMV1oesqU2VTO1Nbko6myWbhpq+kckwST90avjRgYzptYm+YYMwBBk3GQWavZb1NXs2+mm83um1WY
                 bQbyXYEiQWbKZjNM7wJVzjBlYICoPW8CMiXH+LXRR9NwoPkg82xN5IB2MZu2mGabSBnnAbGscYEJj3GV
                 YQAQg/TVGfaA7RI0BsErRjeNeowDgDQPmF9gkmciaJxtArGjzrAKCGWNpYAQAL0kOBE=''')
 
     def test_mixed_linear(self):
-        ndofs, error, lhs = main(nrefine=2, btype='h-std', etype='mixed', degree=1)
+        error, u = main(nrefine=2, btype='h-std', etype='mixed', degree=1)
         with self.subTest('degrees of freedom'):
-            self.assertEqual(ndofs, 34)
+            self.assertEqual(len(u), 34)
         with self.subTest('L2-error'):
             self.assertAlmostEqual(error[0], 0.00450, places=5)
         with self.subTest('H1-error'):
             self.assertAlmostEqual(error[1], 0.11683, places=5)
         with self.subTest('left-hand side'):
-            self.assertAlmostEqual64(lhs, '''
+            self.assertAlmostEqual64(u, '''
                 eNprMT1u6mQyxUTRzMCUAQhazL6b3jNrMYPxp5iA5FtMD+lcMgDxHa4aXzS+6HDV+fKO85cMnC8zMBzS
                 AQDBThbY''')

examples/burgers.py

@@ -37,29 +37,28 @@ def main(nelems: int, ndims: int, btype: str, degree: int, timescale: float, new
     ns = Namespace()
     ns.x = geom
     ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
-    ns.basis = domain.basis(btype, degree=degree)
-    ns.u = function.dotarg('lhs', ns.basis)
+    ns.add_field(('u', 'v'), domain.basis(btype, degree=degree))
     ns.f = '.5 u^2'
     ns.C = 1
     ns.u0 = 'exp(-25 x_i x_i)'
 
-    res = domain.integral('-∇_0(basis_n) f dV' @ ns, degree=5)
-    res += domain.interfaces.integral('-[basis_n] n_0 ({f} - .5 C [u] n_0) dS' @ ns, degree=degree*2)
-    inertia = domain.integral('basis_n u dV' @ ns, degree=5)
+    res = domain.integral('-∇_0(v) f dV' @ ns, degree=5)
+    res += domain.interfaces.integral('-[v] n_0 ({f} - .5 C [u] n_0) dS' @ ns, degree=degree*2)
+    inertia = domain.integral('v u dV' @ ns, degree=5)
 
     sqr = domain.integral('(u - u0)^2 dV' @ ns, degree=5)
-    lhs0 = solver.optimize('lhs', sqr)
+    args0 = solver.optimize('u,', sqr)
 
     timestep = timescale/nelems
     bezier = domain.sample('bezier', 7)
-    with treelog.iter.plain('timestep', solver.impliciteuler('lhs', res, inertia, timestep=timestep, arguments=dict(lhs=lhs0), newtontol=newtontol)) as steps:
-        for itime, lhs in enumerate(steps):
-            x, u = bezier.eval(['x_i', 'u'] @ ns, lhs=lhs)
-            export.triplot('solution.png', x, u, tri=bezier.tri, hull=bezier.hull, clim=(0, 1))
+    with treelog.iter.plain('timestep', solver.impliciteuler('u:v', res, inertia, timestep=timestep, arguments=args0, newtontol=newtontol)) as steps:
+        for itime, args in enumerate(steps):
+            xsmp, usmp = bezier.eval(['x_i', 'u'] @ ns, **args)
+            export.triplot('solution.png', xsmp, usmp, tri=bezier.tri, hull=bezier.hull, clim=(0, 1))
             if itime * timestep >= endtime:
                 break
 
-    return lhs
+    return args['u']
 
 # If the script is executed (as opposed to imported), :func:`nutils.cli.run`
 # calls the main function with arguments provided from the command line. For
@@ -80,38 +79,38 @@ def main(nelems: int, ndims: int, btype: str, degree: int, timescale: float, new
 class test(testing.TestCase):
 
     def test_1d_p0(self):
-        lhs = main(ndims=1, nelems=10, timescale=.1, btype='discont', degree=0, endtime=.01, newtontol=1e-5)
-        self.assertAlmostEqual64(lhs, '''
+        u = main(ndims=1, nelems=10, timescale=.1, btype='discont', degree=0, endtime=.01, newtontol=1e-5)
+        self.assertAlmostEqual64(u, '''
             eNrz1ttqGGOiZSZlrmbuZdZgcsEwUg8AOqwFug==''')
 
     def test_1d_p1(self):
-        lhs = main(ndims=1, nelems=10, timescale=.1, btype='discont', degree=1, endtime=.01, newtontol=1e-5)
-        self.assertAlmostEqual64(lhs, '''
+        u = main(ndims=1, nelems=10, timescale=.1, btype='discont', degree=1, endtime=.01, newtontol=1e-5)
+        self.assertAlmostEqual64(u, '''
             eNrbocann6u3yqjTyMLUwfSw2TWzKPNM8+9mH8wyTMNNZxptMirW49ffpwYAI6cOVA==''')
 
     def test_1d_p2(self):
-        lhs = main(ndims=1, nelems=10, timescale=.1, btype='discont', degree=2, endtime=.01, newtontol=1e-5)
-        self.assertAlmostEqual64(lhs, '''
+        u = main(ndims=1, nelems=10, timescale=.1, btype='discont', degree=2, endtime=.01, newtontol=1e-5)
+        self.assertAlmostEqual64(u, '''
             eNrr0c7SrtWfrD/d4JHRE6Ofxj6mnqaKZofNDpjZmQeYB5pHmL8we23mb5ZvWmjKY/LV6KPRFIMZ+o36
             8dp92gCxZxZG''')
 
     def test_1d_p1_legendre(self):
-        lhs = main(ndims=1, nelems=10, timescale=.1, btype='legendre', degree=1, endtime=.01, newtontol=1e-5)
-        self.assertAlmostEqual64(lhs, '''
+        u = main(ndims=1, nelems=10, timescale=.1, btype='legendre', degree=1, endtime=.01, newtontol=1e-5)
+        self.assertAlmostEqual64(u, '''
             eNrbpbtGt9VQyNDfxMdYzczERNZczdjYnOdsoNmc01kmE870Gj49t0c36BIAAhsO1g==''')
 
     def test_1d_p2_legendre(self):
-        lhs = main(ndims=1, nelems=10, timescale=.1, btype='legendre', degree=2, endtime=.01, newtontol=1e-5)
-        self.assertAlmostEqual64(lhs, '''
+        u = main(ndims=1, nelems=10, timescale=.1, btype='legendre', degree=2, endtime=.01, newtontol=1e-5)
+        self.assertAlmostEqual64(u, '''
             eNoBPADD/8ot2y2/K4UxITFFLk00RTNNLyY2KzTTKx43QjOOzzM3Ss0pz1A2qsvhKGk0jsyXL48xzc5j
             LswtIdLIK5SlF78=''')
 
     def test_2d_p1(self):
-        lhs = main(ndims=2, nelems=4, timescale=.1, btype='discont', degree=1, endtime=.01, newtontol=1e-5)
+        u = main(ndims=2, nelems=4, timescale=.1, btype='discont', degree=1, endtime=.01, newtontol=1e-5)
         import os
         if os.environ.get('NUTILS_TENSORIAL'):
-            lhs = lhs.reshape(4, 2, 4, 2).transpose(0, 2, 1, 3).ravel()
-        self.assertAlmostEqual64(lhs, '''
+            u = u.reshape(4, 2, 4, 2).transpose(0, 2, 1, 3).ravel()
+        self.assertAlmostEqual64(u, '''
             eNoNyKENhEAQRuGEQsCv2SEzyQZHDbRACdsDJNsBjqBxSBxBHIgJ9xsqQJ1Drro1L1/eYBZceGz8njrR
             yacm8UQLBvPYCw1airpyUVYSJLhKijK4IC01WDnqqxvX8OTl427aU73sctPGr3qqceBnRzOjo0xy9JpJ
             R73m6R6YMZo/Q+FCLQ==''')

examples/cahnhilliard.py

@@ -151,22 +151,19 @@ def main(size: Length, epsilon: Length, mobility: Time/Density, stens: Tension,
 
     bezier = domain.sample('bezier', 5)  # sample for surface plots
     grid = domain.locate(geom, numeric.simplex_grid([1, 1], 1/40), maxdist=1/nelems, skip_missing=True, tol=1e-5)  # sample for quivers
-
-    φbasis = ηbasis = domain.basis('std', degree=degree)
-    ηbasis *= stens / epsilon  # basis scaling to give η the required unit
+    basis = domain.basis('std', degree=degree)
 
     ns = Namespace()
     ns.x = size * geom
     ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
+    ns.add_field(('φ', 'φ0'), basis)
+    ns.add_field('η', basis * stens / epsilon) # basis scaling to give η the required unit
     ns.ε = epsilon
     ns.σ = stens
-    ns.φ = function.dotarg('φ', φbasis)
     ns.σmean = (wtensp + wtensn) / 2
     ns.σdiff = (wtensp - wtensn) / 2
     ns.σwall = 'σmean + φ σdiff'
-    ns.φ0 = function.dotarg('φ0', φbasis)
     ns.dφ = 'φ - φ0'
-    ns.η = function.dotarg('η', ηbasis)
     ns.ψ = '.25 (φ^2 - 1)^2'
     ns.δψ = stab.value
     ns.M = mobility
@@ -179,31 +176,31 @@ def main(size: Length, epsilon: Length, mobility: Time/Density, stens: Tension,
     nrg = nrg_mix + nrg_iface + nrg_wall + domain.integral('(δψ σ / ε - η dφ + .5 dt J_k ∇_k(η)) dV' @ ns, degree=7)
 
     numpy.random.seed(seed)
-    state = dict(φ=numpy.random.normal(0, .5, φbasis.shape))  # initial condition
+    args = dict(φ=numpy.random.normal(0, .5, basis.shape)) # initial condition
 
     with log.iter.fraction('timestep', range(round(endtime / timestep))) as steps:
         for istep in steps:
 
-            E = numpy.stack(function.eval([nrg_mix, nrg_iface, nrg_wall], **state))
+            E = numpy.stack(function.eval([nrg_mix, nrg_iface, nrg_wall], **args))
             log.user('energy: {0:,.0μJ/m} ({1[0]:.0f}% mixture, {1[1]:.0f}% interface, {1[2]:.0f}% wall)'.format(numpy.sum(E), 100*E/numpy.sum(E)))
 
-            state['φ0'] = state['φ']
-            state = solver.optimize(['φ', 'η'], nrg / tol, arguments=state, tol=1)
+            args['φ0'] = args['φ']
+            args = solver.optimize(['φ', 'η'], nrg / tol, arguments=args, tol=1)
 
             with export.mplfigure('phase.png') as fig:
                 ax = fig.add_subplot(aspect='equal', xlabel='[mm]', ylabel='[mm]')
-                x, φ = bezier.eval(['x_i', 'φ'] @ ns, **state)
+                x, φ = bezier.eval(['x_i', 'φ'] @ ns, **args)
                 im = ax.tripcolor(*(x/'mm').T, bezier.tri, φ, shading='gouraud', cmap='bwr')
                 im.set_clim(-1, 1)
                 fig.colorbar(im)
                 if showflux:
-                    x, J = grid.eval(['x_i', 'J_i'] @ ns, **state)
+                    x, J = grid.eval(['x_i', 'J_i'] @ ns, **args)
                     log.info('largest flux: {:.1mm/h}'.format(numpy.max(numpy.hypot(J[:, 0], J[:, 1]))))
                     ax.quiver(*(x/'mm').T, *(J/'m/s').T, color='r')
                     ax.quiver(*(x/'mm').T, *-(J/'m/s').T, color='b')
                 ax.autoscale(enable=True, axis='both', tight=True)
 
-    return state
+    return args
 
 # If the script is executed (as opposed to imported), :func:`nutils.cli.run`
 # calls the main function with arguments provided from the command line.
@@ -222,39 +219,39 @@ def main(size: Length, epsilon: Length, mobility: Time/Density, stens: Tension,
 class test(testing.TestCase):
 
     def test_initial(self):
-        state = main(size=parse('10cm'), epsilon=parse('5cm'), mobility=parse('1μL*s/kg'),
+        args = main(size=parse('10cm'), epsilon=parse('5cm'), mobility=parse('1μL*s/kg'),
                      stens=parse('50mN/m'), wtensn=parse('30mN/m'), wtensp=parse('20mN/m'), nelems=3,
                      etype='square', degree=2, timestep=parse('1h'), tol=parse('1nJ/m'),
                      endtime=parse('1h'), seed=0, circle=False, stab=stab.linear, showflux=True)
         with self.subTest('concentration'):
-            self.assertAlmostEqual64(state['φ0'], '''
+            self.assertAlmostEqual64(args['φ0'], '''
                 eNoBYgCd/xM3LjTtNYs3MDcUyt41uc14zjo0LzKzNm812jFhNNMzwDYgzbMzV8o0yCM1rzWeypE3Tcnx
                 L07NzTa4NlMyETREyrPIGMxYMl82VDbjy1/M8clZyf3IRjday6XLmMl6NRnJMF4tqQ==''')
 
     def test_square(self):
-        state = main(size=parse('10cm'), epsilon=parse('5cm'), mobility=parse('1μL*s/kg'),
+        args = main(size=parse('10cm'), epsilon=parse('5cm'), mobility=parse('1μL*s/kg'),
                      stens=parse('50mN/m'), wtensn=parse('30mN/m'), wtensp=parse('20mN/m'), nelems=3,
                      etype='square', degree=2, timestep=parse('1h'), tol=parse('1nJ/m'),
                      endtime=parse('2h'), seed=0, circle=False, stab=stab.linear, showflux=True)
         with self.subTest('concentration'):
-            self.assertAlmostEqual64(state['φ'], '''
+            self.assertAlmostEqual64(args['φ'], '''
                 eNoBYgCd/zE1EjX1NaA2+TXiMxkz0TS9NL01ajaRNZoxYNElNRM1LDUlNZQw0cqgysI1nTWcNN4xLsuk
                 ybDJvDWaNTQ07s7nysnJ6ckPNQY1CzNozKjK58kOysQ0zTQKM73M3coVyjfKR9cuPg==''')
         with self.subTest('chemical-potential'):
-            self.assertAlmostEqual64(state['η'], '''
+            self.assertAlmostEqual64(args['η'], '''
                 eNoBYgCd/3TIkccBNkQ6IDqIN3/MF8cSx9Y02TmdOVHLMcecxxLIEjQUOAHOa8a1xWw3izb2M9UzPMc0
                 xmnGpzibODY34tETyJHHp8hbyWU2xzZTydfIOsrNyo3Gi8jCyyXIm8hkzD3K1IAxtQ==''')
 
     def test_mixedcircle(self):
-        state = main(size=parse('10cm'), epsilon=parse('5cm'), mobility=parse('1μL*s/kg'),
+        args = main(size=parse('10cm'), epsilon=parse('5cm'), mobility=parse('1μL*s/kg'),
                      stens=parse('50mN/m'), wtensn=parse('30mN/m'), wtensp=parse('20mN/m'), nelems=3,
                      etype='mixed', degree=2, timestep=parse('1h'), tol=parse('1nJ/m'),
                      endtime=parse('2h'), seed=0, circle=True, stab=stab.linear, showflux=True)
         with self.subTest('concentration'):
-            self.assertAlmostEqual64(state['φ'], '''
+            self.assertAlmostEqual64(args['φ'], '''
                 eNoBYgCd/w01AjX+NAw1IjXTNMw0iTRPNDI1vDQcNTk0uzJ9NFM0HS4P0SbMcssOy0wzZjNw0b0zljHK
                 z6U0ps8zM/LPjspVypDKUsuLzk3MgM3OzYnN7s/61KfP2zH4MADNhst3z7DMoBcvyQ==''')
         with self.subTest('chemical-potential'):
-            self.assertAlmostEqual64(state['η'], '''
+            self.assertAlmostEqual64(args['η'], '''
                 eNoBYgCd/+s1Bcp4ztI3gjYFyZk4YzVjyfA2AzdAMj032zfLNTE4fMm7yLnGisbqxZPJ2MsfyD81csiv
                 x+E5xDhjOJA3msZ1xZTFa8ddx/fG88eCx73H1MieM/c0WDihMUrLvMYZNpvIrWQ0sw==''')