adaptivity.pyΒΆ

In this script we solve the Laplace problem on a unit square that has the bottom-right quadrant removed (a.k.a. an L-shaped domain) with Dirichlet boundary conditions matching the harmonic function

\[\sqrt[3]{x^2 + y^2} \cos\left(\tfrac23 \arctan\frac{y+x}{y-x}\right),\]

shifted by 0.5 such that the origin coincides with the middle of the unit square. This variation of a well known benchmark problem is known to converge suboptimally under uniform refinement due to a singular gradient in the reentrant corner. This script demonstrates that optimal convergence can be restored by using adaptive refinement.

15from nutils import mesh, function, solver, util, export, cli, testing
16import numpy, treelog
17
18def main(etype:str, btype:str, degree:int, nrefine:int):
19  '''
20  Adaptively refined Laplace problem on an L-shaped domain.
21
22  .. arguments::
23
24     etype [square]
25       Type of elements (square/triangle/mixed).
26     btype [h-std]
27       Type of basis function (h/th-std/spline), with availability depending on
28       the configured element type.
29     degree [2]
30       Polynomial degree
31     nrefine [5]
32       Number of refinement steps to perform.
33  '''
34
35  domain, geom = mesh.unitsquare(2, etype)
36
37  x, y = geom - .5
38  exact = (x**2 + y**2)**(1/3) * function.cos(function.arctan2(y+x, y-x) * (2/3))
39  domain = domain.trim(exact-1e-15, maxrefine=0)
40  linreg = util.linear_regressor()
41
42  with treelog.iter.fraction('level', range(nrefine+1)) as lrange:
43    for irefine in lrange:
44
45      if irefine:
46        refdom = domain.refined
47        ns.refbasis = refdom.basis(btype, degree=degree)
48        indicator = refdom.integral('refbasis_n,k u_,k d:x' @ ns, degree=degree*2).eval(lhs=lhs)
49        indicator -= refdom.boundary.integral('refbasis_n u_,k n_k d:x' @ ns, degree=degree*2).eval(lhs=lhs)
50        supp = ns.refbasis.get_support(indicator**2 > numpy.mean(indicator**2))
51        domain = domain.refined_by(ns.refbasis.transforms[supp])
52
53      ns = function.Namespace()
54      ns.x = geom
55      ns.basis = domain.basis(btype, degree=degree)
56      ns.u = 'basis_n ?lhs_n'
57      ns.du = ns.u - exact
58
59      sqr = domain.boundary['trimmed'].integral('u^2 d:x' @ ns, degree=degree*2)
60      cons = solver.optimize('lhs', sqr, droptol=1e-15)
61
62      sqr = domain.boundary.integral('du^2 d:x' @ ns, degree=7)
63      cons = solver.optimize('lhs', sqr, droptol=1e-15, constrain=cons)
64
65      res = domain.integral('basis_n,k u_,k d:x' @ ns, degree=degree*2)
66      lhs = solver.solve_linear('lhs', res, constrain=cons)
67
68      ndofs = len(ns.basis)
69      error = domain.integral('<du^2, du_,k du_,k>_i d:x' @ ns, degree=7).eval(lhs=lhs)**.5
70      rate, offset = linreg.add(numpy.log(len(ns.basis)), numpy.log(error))
71      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))
72
73      bezier = domain.sample('bezier', 9)
74      x, u, du = bezier.eval(['x_i', 'u', 'du'] @ ns, lhs=lhs)
75      export.triplot('sol.png', x, u, tri=bezier.tri, hull=bezier.hull)
76      export.triplot('err.png', x, du, tri=bezier.tri, hull=bezier.hull)
77
78  return ndofs, error, 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 perform four refinement steps with quadratic basis functions starting from a triangle mesh run python3 adaptivity.py etype=triangle degree=2 nrefine=4 (view log).

86if __name__ == '__main__':
87  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.

 95class test(testing.TestCase):
 96
 97  @testing.requires('matplotlib')
 98  def test_square_quadratic(self):
 99    ndofs, error, lhs = main(nrefine=2, btype='h-std', etype='square', degree=2)
100    with self.subTest('degrees of freedom'):
101      self.assertEqual(ndofs, 149)
102    with self.subTest('L2-error'):
103      self.assertAlmostEqual(error[0], 0.00065, places=5)
104    with self.subTest('H1-error'):
105      self.assertAlmostEqual(error[1], 0.03461, places=5)
106    with self.subTest('left-hand side'): self.assertAlmostEqual64(lhs, '''
107      eNo1j6FrQmEUxT8RBi4KllVfMsl3z/nK4zEmLC6bhsKCw2gSw5IPFsymGbZiWnr+By8Ii7Yhsk3BMtC4
108      Z9sJ223ncs85vzvmM9+Yhix8hDIjtnkdHqQSdDDDj1Qajr5qPXN/07MZ2vI4V7UOIvmdO/oEZY45xYDn
109      oR7ikLHAHVpcs2A1TLhChDO+MOeWt5xjYzm6fOQrGxxiZPeoMGaf37hCyU72hB0u6PglPcQcKxRI/KUd
110      7AYLvMPpsqGkCTPumzWf+qV92kKevjK36ozDP/FSnh1iteWiqWuf+oMaKuyKaC1i52rKPokiF2WLA/20
111      bya+ZCPbWKRPpvgFaedebw==''')
112
113  @testing.requires('matplotlib')
114  def test_triangle_quadratic(self):
115    ndofs, error, lhs = main(nrefine=2, btype='h-std', etype='triangle', degree=2)
116    with self.subTest('degrees of freedom'):
117      self.assertEqual(ndofs, 98)
118    with self.subTest('L2-error'):
119      self.assertAlmostEqual(error[0], 0.00138, places=5)
120    with self.subTest('H1-error'):
121      self.assertAlmostEqual(error[1], 0.05324, places=5)
122    with self.subTest('left-hand side'): self.assertAlmostEqual64(lhs, '''
123      eNprMV1oesqU2VTO1Nbko6myWbhpq+kckwST90avjRgYzptYm+YYMwBBk3GQWavZb1NXs2+mm83um1WY
124      bQbyXYEiQWbKZjNM7wJVzjBlYICoPW8CMiXH+LXRR9NwoPkg82xN5IB2MZu2mGabSBnnAbGscYEJj3GV
125      YQAQg/TVGfaA7RI0BsErRjeNeowDgDQPmF9gkmciaJxtArGjzrAKCGWNpYAQAL0kOBE=''')
126
127  @testing.requires('matplotlib')
128  def test_mixed_linear(self):
129    ndofs, error, lhs = main(nrefine=2, btype='h-std', etype='mixed', degree=1)
130    with self.subTest('degrees of freedom'):
131      self.assertEqual(ndofs, 34)
132    with self.subTest('L2-error'):
133      self.assertAlmostEqual(error[0], 0.00450, places=5)
134    with self.subTest('H1-error'):
135      self.assertAlmostEqual(error[1], 0.11683, places=5)
136    with self.subTest('left-hand side'): self.assertAlmostEqual64(lhs, '''
137      eNprMT1u6mQyxUTRzMCUAQhazL6b3jNrMYPxp5iA5FtMD+lcMgDxHa4aXzS+6HDV+fKO85cMnC8zMBzS
138      AQDBThbY''')