From 682a46767f2d039f28b34356ffe58c5f45b74baa Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 18 Jul 2024 10:46:56 +0100 Subject: [PATCH] Test high-order Argyris and HCT (#3658) * Test high-order Argyris and HCT --- tests/regression/test_helmholtz_zany.py | 40 +++++++------- tests/regression/test_projection_zany.py | 67 +++++++++++++++++++----- 2 files changed, 76 insertions(+), 31 deletions(-) diff --git a/tests/regression/test_helmholtz_zany.py b/tests/regression/test_helmholtz_zany.py index 1325fa31a4..7a5d7fffe0 100644 --- a/tests/regression/test_helmholtz_zany.py +++ b/tests/regression/test_helmholtz_zany.py @@ -20,7 +20,7 @@ def helmholtz(n, el_type, degree, perturb): mesh = UnitSquareMesh(2**n, 2**n) if perturb: - V = FunctionSpace(mesh, mesh.coordinates.ufl_element()) + V = mesh.coordinates.function_space() eps = Constant(1 / 2**(n+1)) x, y = SpatialCoordinate(mesh) @@ -30,38 +30,40 @@ def helmholtz(n, el_type, degree, perturb): V = FunctionSpace(mesh, el_type, degree) x = SpatialCoordinate(V.mesh()) + degree = V.ufl_element().degree() # Define variational problem - lmbda = 1 + k = Constant(3) + u_exact = cos(x[0]*pi*k) * cos(x[1]*pi*k) + + lmbda = Constant(1) u = TrialFunction(V) v = TestFunction(V) - f = Function(V) - f.project((1+8*pi*pi)*cos(x[0]*pi*2)*cos(x[1]*pi*2)) - a = (inner(grad(u), grad(v)) + lmbda * inner(u, v)) * dx - L = inner(f, v) * dx + a = inner(grad(u), grad(v))*dx + lmbda * inner(u, v)*dx + L = a(v, u_exact) # Compute solution - assemble(a) - assemble(L) sol = Function(V) - solve(a == L, sol, solver_parameters={'ksp_type': 'preonly', 'pc_type': 'lu'}) + solve(a == L, sol) - # Analytical solution - f.project(cos(x[0]*pi*2)*cos(x[1]*pi*2)) - return sqrt(assemble(inner(sol - f, sol - f) * dx)) + # Return the H1 norm + return sqrt(assemble(a(sol - u_exact, sol - u_exact))) -# Test convergence on Hermite, Bell, and Argyris +# Test H1 convergence on Hermite, HCT, Bell, and Argyris # Morley is omitted since it only can be used on 4th-order problems. # It is, somewhat oddly, a suitable C^1 nonconforming element but # not a suitable C^0 nonconforming one. @pytest.mark.parametrize(('el', 'deg', 'convrate'), - [('Hermite', 3, 3.8), - ('HCT', 3, 3.8), - ('Bell', 5, 4.8), - ('Argyris', 5, 4.8)]) + [('Hermite', 3, 3), + ('HCT', 3, 3), + ('HCT', 4, 4), + ('Bell', 5, 4), + ('Argyris', 5, 5), + ('Argyris', 6, 6)]) @pytest.mark.parametrize("perturb", [False, True], ids=["Regular", "Perturbed"]) def test_firedrake_helmholtz_scalar_convergence(el, deg, convrate, perturb): - diff = np.array([helmholtz(i, el, deg, perturb) for i in range(1, 4)]) + l = 4 + diff = np.array([helmholtz(i, el, deg, perturb) for i in range(l, l+2)]) conv = np.log2(diff[:-1] / diff[1:]) - assert (np.array(conv) > convrate).all() + assert (np.array(conv) > convrate - 0.3).all() diff --git a/tests/regression/test_projection_zany.py b/tests/regression/test_projection_zany.py index 486dac72c0..db94c9dd46 100644 --- a/tests/regression/test_projection_zany.py +++ b/tests/regression/test_projection_zany.py @@ -2,7 +2,7 @@ on the unit square of a function - f(x, y) = (1.0 + 8.0*pi**2)*cos(x[0]*2*pi)*cos(x[1]*2*pi) + f(x, y) = LegendreP(degree + 2, x + y) using elements with nonstandard pullbacks """ @@ -13,10 +13,22 @@ from firedrake import * +relative_magnitudes = lambda x: np.array(x)[1:] / np.array(x)[:-1] +convergence_orders = lambda x: -np.log2(relative_magnitudes(x)) + + +def jrc(a, b, n): + """Jacobi recurrence coefficients""" + an = (2*n+1+a+b)*(2*n+2+a+b) / (2*(n+1)*(n+1+a+b)) + bn = (a+b)*(a-b)*(2*n+1+a+b) / (2*(n+1)*(n+1+a+b)*(2*n+a+b)) + cn = (n+a)*(n+b)*(2*n+2+a+b) / ((n+1)*(n+1+a+b)*(2*n+a+b)) + return an, bn, cn + + @pytest.fixture def hierarchy(request): msh = UnitSquareMesh(2, 2) - return MeshHierarchy(msh, 2) + return MeshHierarchy(msh, 4) def do_projection(mesh, el_type, degree): @@ -27,16 +39,23 @@ def do_projection(mesh, el_type, degree): v = TestFunction(V) x = SpatialCoordinate(mesh) - f = np.prod([sin((1+i) * x[i]*pi) for i in range(len(x))]) - f = f * Constant(np.ones(u.ufl_shape)) - + m = degree + 2 + z = sum(x) + p0 = 1 + p1 = z + for n in range(1, m): + an, bn, cn = jrc(0, 0, n) + ptemp = p1 + p1 = (an * z + bn) * p1 - cn * p0 + p0 = ptemp + + f = p1 * Constant(np.ones(u.ufl_shape)) a = inner(u, v) * dx L = inner(f, v) * dx # Compute solution fh = Function(V) - solve(a == L, fh, solver_parameters={'ksp_type': 'preonly', 'pc_type': 'lu'}) - + solve(a == L, fh) return sqrt(assemble(inner(fh - f, fh - f) * dx)) @@ -44,11 +63,35 @@ def do_projection(mesh, el_type, degree): [('Johnson-Mercier', 1, 1.8), ('Morley', 2, 2.4), ('HCT-red', 3, 2.7), - ('HCT', 3, 3), - ('Hermite', 3, 3), - ('Bell', 5, 4), - ('Argyris', 5, 4.9)]) + ('HCT', 3, 3.7), + ('HCT', 4, 4.8), + ('Hermite', 3, 3.8), + ('Bell', 5, 4.7), + ('Argyris', 5, 5.8), + ('Argyris', 6, 6.7)]) def test_projection_zany_convergence_2d(hierarchy, el, deg, convrate): - diff = np.array([do_projection(m, el, deg) for m in hierarchy]) + diff = np.array([do_projection(m, el, deg) for m in hierarchy[2:]]) conv = np.log2(diff[:-1] / diff[1:]) assert (np.array(conv) > convrate).all() + + +@pytest.mark.parametrize(('element', 'degree'), + [('HCT-red', 3), + ('HCT', 3), + ('HCT', 4), + ('Argyris', 5), + ('Argyris', 6)]) +def test_mass_conditioning(element, degree, hierarchy): + mass_cond = [] + for msh in hierarchy[1:4]: + V = FunctionSpace(msh, element, degree) + u = TrialFunction(V) + v = TestFunction(V) + a = inner(u, v)*dx + B = assemble(a, mat_type="aij").M.handle + A = B.convert("dense").getDenseArray() + kappa = np.linalg.cond(A) + + mass_cond.append(kappa) + + assert max(relative_magnitudes(mass_cond)) < 1.1