Finite element exterior calculus (FEEC) implementation - problems.

[ethrelfall]

Dear community,

I'm trying to implement some of the methods in D.N. Arnold's textbook on FEEC, in Firedrake.

(I think) I have set up a 2D Hodge Laplacian problem for a vector field, but the type of function spaces Arnold recommends do not seem to work well - the vector field output looks "messy" (and I cannot get higher polynomial orders than one to converge) and just using ordinary Lagrange elements appears to work better.

Can anyone say anything to this, please, e.g. are there any obvious mistakes (code, or maths) in my file (I only started using Firedrake a week ago, so there probably are).

Thanks,

Ed.

Firedrake script, which I've attempted to comment suitably:

https://github.com/ethrelfall/Finite-element-exterior-calculus/blob/443ad042a461be5d601dfac62c722c92973bada6/FEEC_test_1.py

[wence-]

So the problem is you're trying to solve a 1-form Hodge-Laplacian problem in 2D, the mixed formulation (4.32) of Arnold is:

Find $(\sigma, u) \in H^1 \times H(\text{curl})$ such that

$$ \begin{align} (\sigma, \tau) - (u, \operatorname{grad} \tau) &= 0 \\ (\operatorname{grad} \sigma, v) + (\operatorname{curl} u, \operatorname{curl} v) + (p, v) &= (f, v) \end{align} $$

For all $(\tau, v) \in H^1 \times H(\text{curl})$.

Eq (7.10) of Arnold gives FEEC names for pairs of spaces that are stable for discretisation of the k-form Hodge Laplacian, in your code you wanted to use the $P^-$ spaces: $V^{k-1} := P_r^{-}\Lambda^{k-1} \subset H^1$, $V^k := P_{r}^{-}\Lambda^{k} \subset H(\text{curl})$. Looking at the femtable we can see that in 2D the $V^0$ space is piecewise polynomials (Lagrange elements) of degree $r$, and the $V^1$ space is RT elements of degree $r$.

But we have to be slightly careful, there are two versions of the "RT" elements in 2D, related by a rotation of the dofs by 90 degrees. You want the $H(\text{curl})$ version of the RT elements, which are edge degrees of freedom (tangential components) rather than the "classic" $H(\text{div})$ versions (normal components).

You code asks for the "face" version of the RT elements FunctionSpace(mesh, "RT", 1), rather than the "edge" version FunctionSpace(mesh, "RTE", 1), which was your only mistake.

So this works:

from firedrake import *
mesh = ...
# The V^1 space
SRT = FunctionSpace(mesh, "RTE", 1)
# The V^0 space
SP = FunctionSpace(mesh, "CG", 1)
V = SRT * SP

u, sigma = TrialFunctions(V)
v, tau = TestFunctions(V)

a = (
    inner(sigma, tau)
    - inner(u, grad(tau))
    + inner(grad(sigma), v)
    + inner(curl(u), curl(v))
) * dx

x, y = SpatialCoordinate(mesh)

f = as_vector([-1, 0])
L = inner(f, v) * dx

g = Function(V)
solve(a == L, g)

[ethrelfall]

Thanks so much for this prompt reply. I'm sufficiently new to FEEC that I didn't know there were two forms of the 2D RT element.

I have tried re-running with the V^1 space set to RTE as suggested - the convergence seems very slow. May I ask additionally whether there is a recommendation for the solver parameters?

Thanks again.

[wence-]

For 2D problems if you specify no solver parameters at all you get sparse LU, which will be fast up to a millions dofs or so. Beyond that it depends a bit on what your mesh is. The canonical preconditioners for problems like this are block diagonal inverses of the Riesz map. If you have a geometry which can be captured by a coarse mesh then you can use a mesh hierarchy and geometric multigrid with Arnold Falk Winther overlapping additive Schwarz smoothers for the curl-curl problem. If your geometry isn't amenable to this then you can use (only at lowest order right now), I think, the Hiptmair-Xu style auxiliary space preconditioners available in hypre.

Or you might be able to do monolithic multigrid (again needing a geometric mesh hierarchy).

See the available options that are not just "black box" here: https://www.firedrakeproject.org/preconditioning.html with some more discussion on preconditioners for saddle point systems here: https://www.firedrakeproject.org/demos/saddle_point_systems.py.html

[ethrelfall]

I was running unpreconditioned gmres (which seemed to work for the order-1 "wrong" RT elements) and not getting convergence with RTE but now I see that with the default solver parameters it works just fine. Thanks again for the excellent support.

[colinjcotter]

Also, Ed, are you the same Ed Threlfall who went to Farlingaye High School in the 1990s?

[ethrelfall]

Colin, I try not to think too much about my years at FHS, but yes, it's me. Good to hear from you :)

[colinjcotter]

Wow! Drop me a line by email if you want to catch up. On 27 May 2022, at 21:14, ethrelfall ***@***.***> wrote:  Colin, I try not to think too much about my years at FHS, but yes, it's me. Good to hear from you :)