Code generation error with serendipity elements and a biharmonic term

[pefarrell]

We are trying to implement the C0 interior penalty method for biharmonic problems described in https://doi.org/10.1016/S0045-7825(02)00286-4. It converges very nicely on quads with CG elements.

We would like to try serendipity elements as we expect them to be cheaper. However, something goes wrong in the code generation stack with the cellwise biharmonic term. The following snippet reproduces the problem:

from firedrake import *

mesh = UnitSquareMesh(2, 2, quadrilateral=True)
U = FunctionSpace(mesh, "S", 2)
u = Function(U)
v = TestFunction(U)

E = (
    - 1/2 * u**2 * dx
    + inner(grad(grad(u)), grad(grad(u))) * dx
    )

F = derivative(E, u, v)

bcdata = Constant(1)
bcs = [DirichletBC(U, project(bcdata, U), "on_boundary")]

params = {
        "mat_type": "aij",
        "ksp_type": "preonly",
        "pc_type": "lu",
        }
nvproblem = NonlinearVariationalProblem(F, u, bcs=bcs)
nvsolver  = NonlinearVariationalSolver(nvproblem, solver_parameters=params)
nvsolver.solve()

This yields the backtrace

Traceback (most recent call last):
  File "seren.py", line 19, in <module>
    nvsolver  = NonlinearVariationalSolver(nvproblem, solver_parameters=params)
  File "/scratch/farrellp/install-scripts/firedrake/firedrake-dev-20200608-mpich/src/pyadjoint/pyadjoint/tape.py", line 46, in wrapper
    return function(*args, **kwargs)
  File "/scratch/farrellp/install-scripts/firedrake/firedrake-dev-20200608-mpich/src/firedrake/firedrake/adjoint/variational_solver.py", line 27, in wrapper
    init(self, problem, *args, **kwargs)
  File "/scratch/farrellp/install-scripts/firedrake/firedrake-dev-20200608-mpich/src/firedrake/firedrake/variational_solver.py", line 193, in __init__
    options_prefix=self.options_prefix)
  File "/scratch/farrellp/install-scripts/firedrake/firedrake-dev-20200608-mpich/src/firedrake/firedrake/solving_utils.py", line 140, in __init__
    form_compiler_parameters=self.fcp)
  File "/scratch/farrellp/install-scripts/firedrake/firedrake-dev-20200608-mpich/src/firedrake/firedrake/assemble.py", line 172, in create_assembly_callable
    loops = tuple(loops)
  File "/scratch/farrellp/install-scripts/firedrake/firedrake-dev-20200608-mpich/src/firedrake/firedrake/assemble.py", line 741, in _assemble
    form_compiler_parameters=form_compiler_parameters)
  File "/scratch/farrellp/install-scripts/firedrake/firedrake-dev-20200608-mpich/src/firedrake/firedrake/assemble.py", line 514, in create_parloops
    kernels = tsfc_interface.compile_form(expr, "form", parameters=form_compiler_parameters, diagonal=diagonal)
  File "/scratch/farrellp/install-scripts/firedrake/firedrake-dev-20200608-mpich/src/firedrake/firedrake/tsfc_interface.py", line 218, in compile_form
    number_map, interface, coffee, diagonal).kernels
  File "/scratch/farrellp/install-scripts/firedrake/firedrake-dev-20200608-mpich/src/PyOP2/pyop2/caching.py", line 199, in __new__
    obj = make_obj()
  File "/scratch/farrellp/install-scripts/firedrake/firedrake-dev-20200608-mpich/src/PyOP2/pyop2/caching.py", line 189, in make_obj
    obj.__init__(*args, **kwargs)
  File "/scratch/farrellp/install-scripts/firedrake/firedrake-dev-20200608-mpich/src/firedrake/firedrake/tsfc_interface.py", line 124, in __init__
    tree = tsfc_compile_form(form, prefix=name, parameters=parameters, interface=interface, coffee=coffee, diagonal=diagonal)
  File "/scratch/farrellp/install-scripts/firedrake/firedrake-dev-20200608-mpich/src/tsfc/tsfc/driver.py", line 60, in compile_form
    kernel = compile_integral(integral_data, fd, prefix, parameters, interface=interface, coffee=coffee, diagonal=diagonal)
  File "/scratch/farrellp/install-scripts/firedrake/firedrake-dev-20200608-mpich/src/tsfc/tsfc/driver.py", line 207, in compile_integral
    **config)
  File "/scratch/farrellp/install-scripts/firedrake/firedrake-dev-20200608-mpich/src/tsfc/tsfc/fem.py", line 645, in compile_ufl
    result = map_expr_dags(context.translator, expressions)
  File "/scratch/farrellp/install-scripts/firedrake/firedrake-dev-20200608-mpich/src/ufl/ufl/corealg/map_dag.py", line 73, in map_expr_dags
    r = handlers[v._ufl_typecode_](v)
  File "/scratch/farrellp/install-scripts/firedrake/firedrake-dev-20200608-mpich/src/tsfc/tsfc/ufl_utils.py", line 145, in _modified_terminal
    return self.modified_terminal(o)
  File "/scratch/farrellp/install-scripts/firedrake/firedrake-dev-20200608-mpich/src/tsfc/tsfc/fem.py", line 303, in modified_terminal
    return translate(mt.terminal, mt, self.context)
  File "/usr/lib/python3.6/functools.py", line 807, in wrapper
    return dispatch(args[0].__class__)(*args, **kw)
  File "/scratch/farrellp/install-scripts/firedrake/firedrake-dev-20200608-mpich/src/tsfc/tsfc/fem.py", line 581, in translate_coefficient
    finat_dict = ctx.basis_evaluation(element, mt, entity_id)
  File "/scratch/farrellp/install-scripts/firedrake/firedrake-dev-20200608-mpich/src/tsfc/tsfc/fem.py", line 228, in basis_evaluation
    coordinate_mapping=CoordinateMapping(mt, self))
  File "/scratch/farrellp/install-scripts/firedrake/firedrake-dev-20200608-mpich/src/FInAT/finat/fiat_elements.py", line 78, in basis_evaluation
    space_dimension, value_size, len(ps.points)
ValueError: cannot reshape array of size 8 into shape (8,1,16)

The code compiles fine if one replaces "S" with "CG".

[miklos1]

FIAT.Serendipity is broken for higher derivatives. Minimal failing example:

from pprint import pprint

from FIAT import ufc_cell, make_quadrature, Serendipity

cell = ufc_cell("quadrilateral")
elem = Serendipity(cell, 3)
rule = make_quadrature(cell, 2)
tabs = elem.tabulate(3, rule.get_points())

pprint({alpha: table.shape for alpha, table in tabs.items()})

prints

{(0, 0): (12, 4),
 (0, 1): (12, 4),
 (0, 2): (12,),
 (0, 3): (12,),
 (1, 0): (12, 4),
 (1, 1): (12,),
 (1, 2): (12,),
 (2, 0): (12,),
 (2, 1): (12,),
 (3, 0): (12,)}

The wrong shape then triggers an assertion failure in the FInAT element wrapper of the FIAT element.

[wence-]

Probably this:

diff --git a/FIAT/serendipity.py b/FIAT/serendipity.py
index 3792f70..1f1dbff 100644
--- a/FIAT/serendipity.py
+++ b/FIAT/serendipity.py
@@ -150,6 +150,8 @@ class Serendipity(FiniteElement):
             raise NotImplementedError('no tabulate method for serendipity elements of dimension 1 or less.')
         if dim >= 4:
             raise NotImplementedError('tabulate does not support higher dimensions than 3.')
+        points = np.asarray(points)
+        npoints, pointdim = points.shape
         for o in range(order + 1):
             alphas = mis(dim, o)
             for alpha in alphas:
@@ -160,8 +162,10 @@ class Serendipity(FiniteElement):
                     callable = lambdify(variables[:dim], polynomials, modules="numpy", dummify=True)
                     self.basis[alpha] = polynomials
                     self.basis_callable[alpha] = callable
-                points = np.asarray(points)
-                T = np.asarray(callable(*(points[:, i] for i in range(points.shape[1]))))
+                tabulation = tuple(callable(*(points[:, i]
+                                              for i in range(pointdim))))
+                T = np.asarray([np.broadcast_to(tab, (npoints, ))
+                                for tab in tabulation])
                 phivals[alpha] = T
         return phivals

[wence-]

I think that branch fixes it. @pefarrell can you confirm?

[pefarrell]

It certainly makes the MRE snippet run. I'll ask my student to run our MMS with this patch applied and get back to you.

[pefarrell]

Yep, an MMS for the biharmonic shows the expected convergence rates:

from firedrake import *
import numpy

sp = {"snes_type": "newtonls",
      "snes_convergence_test": "skip",
      "snes_max_it": 2,
      "snes_monitor": None,
      "ksp_type": "gmres",
      "ksp_monitor_true_residual": None,
      "pc_type": "lu",
      "pc_factor_mat_solver_type": "mumps",
      "mat_mumps_icntl_14": 200,}

def error(N):
    mesh = UnitSquareMesh(N, N, quadrilateral=True)
    V = FunctionSpace(mesh, "S", 2)
    u = Function(V)
    v = TestFunction(V)

    X = SpatialCoordinate(mesh)
    k = 1
    u_ex = sin(2*pi*k*X[0]) * cos(2*pi*k*X[1])

    h = avg(CellDiameter(mesh))
    n = FacetNormal(mesh)

    f = div(div(grad(grad(u_ex)))) + u_ex
    g = dot(grad(grad(u_ex)), n)
    F = inner(grad(grad(u)), grad(grad(v)))*dx + inner(u, v)*dx - inner(f, v)*dx - inner(g, grad(v))*ds + h**(-3) * inner(jump(grad(u), n), jump(grad(v), n))*dS

    bc = DirichletBC(V, project(u_ex, V, solver_parameters=sp), "on_boundary")

    solve(F == 0, u, bc, solver_parameters=sp)

    err = errornorm(u_ex, u, "L2")
    return err

errors = []
for N in [10, 20, 40, 80]:
    errors.append(error(N))

print(errors)
def convergence_order(errors):
    from math import log
    out = []
    for (preverror, nexterror) in zip(errors, errors[1:]):
        conv = log(preverror/nexterror, 2)
        out.append(conv)
    return out

conv = convergence_order(errors)
print(conv)

yields

  Residual norms for firedrake_0_ solve.
  0 KSP preconditioned resid norm 5.675064933304e+00 true resid norm 4.584728486039e-02 ||r(i)||/||b|| 1.000000000000e+00
  1 KSP preconditioned resid norm 5.605996400087e-14 true resid norm 8.675951275143e-18 ||r(i)||/||b|| 1.892358795414e-16
  0 SNES Function norm 2.926061223240e+04 
    Residual norms for firedrake_1_ solve.
    0 KSP preconditioned resid norm 4.637503020943e+00 true resid norm 2.926061223240e+04 ||r(i)||/||b|| 1.000000000000e+00
    1 KSP preconditioned resid norm 4.151192971108e-13 true resid norm 2.904764117405e-10 ||r(i)||/||b|| 9.927215788698e-15
  1 SNES Function norm 2.881993078224e-10 
    Residual norms for firedrake_1_ solve.
    0 KSP preconditioned resid norm 1.171427695574e-12 true resid norm 2.881993078224e-10 ||r(i)||/||b|| 1.000000000000e+00
    1 KSP preconditioned resid norm 1.545610282161e-25 true resid norm 3.543360851879e-23 ||r(i)||/||b|| 1.229482776573e-13
  2 SNES Function norm 8.019843441692e-12 
  Residual norms for firedrake_2_ solve.
  0 KSP preconditioned resid norm 1.051317262660e+01 true resid norm 2.465730194105e-02 ||r(i)||/||b|| 1.000000000000e+00
  1 KSP preconditioned resid norm 1.071956664922e-13 true resid norm 4.769524850526e-18 ||r(i)||/||b|| 1.934325524313e-16
  0 SNES Function norm 6.708789181521e+05 
    Residual norms for firedrake_3_ solve.
    0 KSP preconditioned resid norm 9.498624180717e+00 true resid norm 6.708789181521e+05 ||r(i)||/||b|| 1.000000000000e+00
    1 KSP preconditioned resid norm 1.856059729254e-10 true resid norm 1.196926001121e-07 ||r(i)||/||b|| 1.784116281992e-13
  1 SNES Function norm 1.197395028405e-07 
    Residual norms for firedrake_3_ solve.
    0 KSP preconditioned resid norm 3.065653044677e-11 true resid norm 1.197395028405e-07 ||r(i)||/||b|| 1.000000000000e+00
    1 KSP preconditioned resid norm 7.249854561276e-23 true resid norm 3.670434163096e-19 ||r(i)||/||b|| 3.065349426067e-12
  2 SNES Function norm 2.649419675719e-10 
  Residual norms for firedrake_4_ solve.
  0 KSP preconditioned resid norm 2.049704577867e+01 true resid norm 1.263458196617e-02 ||r(i)||/||b|| 1.000000000000e+00
  1 KSP preconditioned resid norm 2.388712764073e-13 true resid norm 1.709292562358e-17 ||r(i)||/||b|| 1.352868315655e-15
  0 SNES Function norm 1.526424092323e+07 
    Residual norms for firedrake_5_ solve.
    0 KSP preconditioned resid norm 1.949092605650e+01 true resid norm 1.526424092323e+07 ||r(i)||/||b|| 1.000000000000e+00
    1 KSP preconditioned resid norm 1.901571392502e-09 true resid norm 2.212415734171e-04 ||r(i)||/||b|| 1.449410910964e-11
  1 SNES Function norm 2.212428144891e-04 
    Residual norms for firedrake_5_ solve.
    0 KSP preconditioned resid norm 1.693401381091e-08 true resid norm 2.212428144891e-04 ||r(i)||/||b|| 1.000000000000e+00
    1 KSP preconditioned resid norm 6.749447836495e-19 true resid norm 9.465038468891e-15 ||r(i)||/||b|| 4.278122428857e-11
  2 SNES Function norm 8.449415673457e-09 
  Residual norms for firedrake_6_ solve.
  0 KSP preconditioned resid norm 4.049730680242e+01 true resid norm 6.375658034144e-03 ||r(i)||/||b|| 1.000000000000e+00
  1 KSP preconditioned resid norm 5.290594861608e-13 true resid norm 3.714227308483e-17 ||r(i)||/||b|| 5.825637586884e-15
  0 SNES Function norm 3.459043547658e+08 
    Residual norms for firedrake_7_ solve.
    0 KSP preconditioned resid norm 3.949449791920e+01 true resid norm 3.459043547658e+08 ||r(i)||/||b|| 1.000000000000e+00
    1 KSP preconditioned resid norm 3.408828072800e-07 true resid norm 4.264925050868e-01 ||r(i)||/||b|| 1.232978131702e-09
  1 SNES Function norm 4.264924800584e-01 
    Residual norms for firedrake_7_ solve.
    0 KSP preconditioned resid norm 4.379326409673e-07 true resid norm 4.264924800584e-01 ||r(i)||/||b|| 1.000000000000e+00
    1 KSP preconditioned resid norm 8.130152694047e-16 true resid norm 1.972246078777e-09 ||r(i)||/||b|| 4.624339633155e-09
  2 SNES Function norm 2.836560353563e-07 
[0.003333349959638379, 0.0006761160289255834, 0.0001526720275078284, 3.6493096778949144e-05]
[2.301630034996994, 2.146835093558506, 2.0647402691534418]

Interesting observation: I can get higher-order convergence with CGk if I use h^{-k-1} in the penalty term, but that doesn't work for serendipity elements. I would expect this is something to be analysed rather than a bug.

[wence-]

Interesting observation: I can get higher-order convergence with CGk if I use h^{-k-1} in the penalty term, but that doesn't work for serendipity elements. I would expect this is something to be analysed rather than a bug.

I remember the convergence of the IP scheme for biharmonic that rob and I did in the zany element paper was quite sensitive to the weighting of the penalty parameter: it might be that you're losing coercivity/continuity? (try with alpha h^{-k-1} with alpha > 1?)

[pefarrell]

Great point! You're right, but with the direction wrong. Setting alpha < 1 recovers kth order convergence:

from firedrake import *
import numpy

sp = {"snes_type": "newtonls",
      "snes_convergence_test": "skip",
      "snes_max_it": 2,
      "snes_monitor": None,
      "ksp_type": "gmres",
      "ksp_monitor_true_residual": None,
      "pc_type": "lu",
      "pc_factor_mat_solver_type": "mumps",
      "mat_mumps_icntl_14": 200,}

def error(N):
    mesh = UnitSquareMesh(N, N, quadrilateral=True)
    degree = 4
    V = FunctionSpace(mesh, "S", degree)
    u = Function(V)
    v = TestFunction(V)

    X = SpatialCoordinate(mesh)
    k = 1
    u_ex = sin(2*pi*k*X[0]) * cos(2*pi*k*X[1])

    h = avg(CellDiameter(mesh))
    alpha = Constant(1/100)
    n = FacetNormal(mesh)

    f = div(div(grad(grad(u_ex)))) + u_ex
    g = dot(grad(grad(u_ex)), n)
    F = inner(grad(grad(u)), grad(grad(v)))*dx + inner(u, v)*dx - inner(f, v)*dx - inner(g, grad(v))*ds + alpha * h**(-(degree+1)) * inner(jump(grad(u), n), jump(grad(v), n))*dS

    bc = DirichletBC(V, project(u_ex, V, solver_parameters=sp), "on_boundary")

    solve(F == 0, u, bc, solver_parameters=sp)

    err = errornorm(u_ex, u, "L2")
    return err

errors = []
for N in [10, 20, 40, 80]:
    errors.append(error(N))

print(errors)
def convergence_order(errors):
    from math import log
    out = []
    for (preverror, nexterror) in zip(errors, errors[1:]):
        conv = log(preverror/nexterror, 2)
        out.append(conv)
    return out

conv = convergence_order(errors)
print(conv)