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
16from nutils.expression_v2 import Namespace
17import numpy
18import treelog
19
20
21def main(etype: str, btype: str, degree: int, nrefine: int):
22 '''
23 Adaptively refined Laplace problem on an L-shaped domain.
24
25 .. arguments::
26
27 etype [square]
28 Type of elements (square/triangle/mixed).
29 btype [h-std]
30 Type of basis function (h/th-std/spline), with availability depending on
31 the configured element type.
32 degree [2]
33 Polynomial degree
34 nrefine [5]
35 Number of refinement steps to perform.
36 '''
37
38 domain, geom = mesh.unitsquare(2, etype)
39
40 x, y = geom - .5
41 exact = (x**2 + y**2)**(1/3) * function.cos(function.arctan2(y+x, y-x) * (2/3))
42 domain = domain.trim(exact-1e-15, maxrefine=0)
43 linreg = util.linear_regressor()
44
45 with treelog.iter.fraction('level', range(nrefine+1)) as lrange:
46 for irefine in lrange:
47
48 if irefine:
49 refdom = domain.refined
50 ns.refbasis = refdom.basis(btype, degree=degree)
51 indicator = refdom.integral('∇_k(refbasis_n) ∇_k(u) dV' @ ns, degree=degree*2).eval(lhs=lhs)
52 indicator -= refdom.boundary.integral('refbasis_n ∇_k(u) n_k dS' @ ns, degree=degree*2).eval(lhs=lhs)
53 supp = ns.refbasis.get_support(indicator**2 > numpy.mean(indicator**2))
54 domain = domain.refined_by(refdom.transforms[supp])
55
56 ns = Namespace()
57 ns.x = geom
58 ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
59 ns.basis = domain.basis(btype, degree=degree)
60 ns.u = function.dotarg('lhs', ns.basis)
61 ns.du = ns.u - exact
62
63 sqr = domain.boundary['trimmed'].integral('u^2 dS' @ ns, degree=degree*2)
64 cons = solver.optimize('lhs', sqr, droptol=1e-15)
65
66 sqr = domain.boundary.integral('du^2 dS' @ ns, degree=7)
67 cons = solver.optimize('lhs', sqr, droptol=1e-15, constrain=cons)
68
69 res = domain.integral('∇_k(basis_n) ∇_k(u) dV' @ ns, degree=degree*2)
70 lhs = solver.solve_linear('lhs', res, constrain=cons)
71
72 ndofs = len(ns.basis)
73 error = function.sqrt(domain.integral(['du du dV', '∇_k(du) ∇_k(du) dV'] @ ns, degree=7)).eval(lhs=lhs)
74 rate, offset = linreg.add(numpy.log(len(ns.basis)), numpy.log(error))
75 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))
76
77 bezier = domain.sample('bezier', 9)
78 x, u, du = bezier.eval(['x_i', 'u', 'du'] @ ns, lhs=lhs)
79 export.triplot('sol.png', x, u, tri=bezier.tri, hull=bezier.hull)
80 export.triplot('err.png', x, du, tri=bezier.tri, hull=bezier.hull)
81
82 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).
91if __name__ == '__main__':
92 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.
101class test(testing.TestCase):
102
103 def test_square_quadratic(self):
104 ndofs, error, lhs = main(nrefine=2, btype='h-std', etype='square', degree=2)
105 with self.subTest('degrees of freedom'):
106 self.assertEqual(ndofs, 149)
107 with self.subTest('L2-error'):
108 self.assertAlmostEqual(error[0], 0.00065, places=5)
109 with self.subTest('H1-error'):
110 self.assertAlmostEqual(error[1], 0.03461, places=5)
111 with self.subTest('left-hand side'):
112 self.assertAlmostEqual64(lhs, '''
113 eNo1j6FrQmEUxT8RBi4KllVfMsl3z/nK4zEmLC6bhsKCw2gSw5IPFsymGbZiWnr+By8Ii7Yhsk3BMtC4
114 Z9sJ223ncs85vzvmM9+Yhix8hDIjtnkdHqQSdDDDj1Qajr5qPXN/07MZ2vI4V7UOIvmdO/oEZY45xYDn
115 oR7ikLHAHVpcs2A1TLhChDO+MOeWt5xjYzm6fOQrGxxiZPeoMGaf37hCyU72hB0u6PglPcQcKxRI/KUd
116 7AYLvMPpsqGkCTPumzWf+qV92kKevjK36ozDP/FSnh1iteWiqWuf+oMaKuyKaC1i52rKPokiF2WLA/20
117 bya+ZCPbWKRPpvgFaedebw==''')
118
119 def test_triangle_quadratic(self):
120 ndofs, error, lhs = main(nrefine=2, btype='h-std', etype='triangle', degree=2)
121 with self.subTest('degrees of freedom'):
122 self.assertEqual(ndofs, 98)
123 with self.subTest('L2-error'):
124 self.assertAlmostEqual(error[0], 0.00138, places=5)
125 with self.subTest('H1-error'):
126 self.assertAlmostEqual(error[1], 0.05324, places=5)
127 with self.subTest('left-hand side'):
128 self.assertAlmostEqual64(lhs, '''
129 eNprMV1oesqU2VTO1Nbko6myWbhpq+kckwST90avjRgYzptYm+YYMwBBk3GQWavZb1NXs2+mm83um1WY
130 bQbyXYEiQWbKZjNM7wJVzjBlYICoPW8CMiXH+LXRR9NwoPkg82xN5IB2MZu2mGabSBnnAbGscYEJj3GV
131 YQAQg/TVGfaA7RI0BsErRjeNeowDgDQPmF9gkmciaJxtArGjzrAKCGWNpYAQAL0kOBE=''')
132
133 def test_mixed_linear(self):
134 ndofs, error, lhs = main(nrefine=2, btype='h-std', etype='mixed', degree=1)
135 with self.subTest('degrees of freedom'):
136 self.assertEqual(ndofs, 34)
137 with self.subTest('L2-error'):
138 self.assertAlmostEqual(error[0], 0.00450, places=5)
139 with self.subTest('H1-error'):
140 self.assertAlmostEqual(error[1], 0.11683, places=5)
141 with self.subTest('left-hand side'):
142 self.assertAlmostEqual64(lhs, '''
143 eNprMT1u6mQyxUTRzMCUAQhazL6b3jNrMYPxp5iA5FtMD+lcMgDxHa4aXzS+6HDV+fKO85cMnC8zMBzS
144 AQDBThbY''')