Hodge decomposition on an annular domain.

[agrubertx]

Hi everyone,

I'm new to Firedrake and attempting to compute the Hodge decomposition of a vector field (i.e., one-form) $f$ (left plot below) on a simple annular domain following system (7.1) in AFW (https://www-users.cse.umn.edu/~arnold//papers/acta.pdf). I am working at the $H^1, H(\mathrm{curl})$ end of the complex, solving the mixed system with a combination of CG and RTE elements by using Firedrake's VectorSpaceBasis functionality to remove the Lagrange multiplier for the harmonic part. However, I am struggling to interpret my results in a way that will allow me to extend things to other meshes and topologies. This is due to a combination of my own ignorance about FEEC and Firedrake; I lifted several aspects of the code (including solver options) from this example. The code I am using is at the end of this message (I couldn't figure out how to attach it).

In my understanding, I have done the following. First, I have projected out the harmonic part of $f$ using the cohomology generator $q$, leaving $\bar{f} = f-P_{\mathfrak{H}}f$ to serve as the problem RHS. This removes the need for the LM $p$ as long as Firedrake knows the nullspace for $u$ should be spanned by $q$. Then, I have solved the resulting linear system for a CG function $\sigma$ whose (weak) gradient generates the exact part of $f$, and an RTE vector field $u$ which can be used to extract the coexact part of $f$. Said differently, I have computed the weak solution to the PDE $d\sigma + \delta du + P_{\mathfrak{H}}f = f$ with minimal norm (c.f., system (7.2)), from which I can compute the coexact part of $f$ by either forming the subtraction $\bar{f}-d\sigma$ or directly computing the weak curl-curl $\delta d u$.

With this said, I am a bit uneasy about the results of my computation. I thought that the components $f_d,f_{\delta},f_h$ of my decomposition (the right three plots below) were supposed to be discrete-orthogonal, but I only get that the inner products $(f_d,f_h), (f_\delta,f_h)$ involving the harmonic part are on the order $\sim 10^{-5}$ while $(f_d,f_\delta)$ is on the order $\sim 10^{-8}$. Is this expected? I also see that $(u,q)$ is on the order $\sim 10^{-3}$, but this orthogonality should only hold weakly, so I am less concerned about this.

Beyond these questions, I have a couple of other related things that I am curious about:

1. Is it possible to implement the Hodge decomposition in Firedrake analogously to (7.1)? If so, is there a reason to prefer one formulation over the other?
2. To extend this to other meshes with more complicated cohomology, is it enough to build a VectorSpaceBasis object from discrete generators which are numerically determined through, e.g., solving the Hodge Laplacian eigenvalue problem?

Thanks in advance for your help,
Anthony

from firedrake import *
import matplotlib.pyplot as plt


######## Creating the mesh #########

import gmsh

# Initialize gmsh
gmsh.initialize()

# Create an annular mesh using gmsh
inner_radius = 0.25
outer_radius = 1.0
mesh_size = 0.04
gmsh.model.add("annulus")

# Create the outer and inner disks
outer_disk = gmsh.model.occ.addDisk(0, 0, 0, outer_radius, outer_radius)
inner_disk = gmsh.model.occ.addDisk(0, 0, 0, inner_radius, inner_radius)

# Cut the inner disk from the outer disk to create the annulus
annulus = gmsh.model.occ.cut([(2, outer_disk)], [(2, inner_disk)], removeTool=True)
gmsh.model.occ.synchronize()

# Set mesh size
gmsh.model.mesh.setSize(gmsh.model.getEntities(0), mesh_size)

# Generate the mesh
gmsh.model.mesh.generate(2)
gmsh.write("annulus.msh")

# Finalize gmsh
gmsh.finalize()


######## Solving the PDE #########

# Load Firedrake mesh
mesh = Mesh("annulus.msh")

# The V^1 space
V1 = FunctionSpace(mesh, "RTE", 3)

# The V^0 space
V0 = FunctionSpace(mesh, "CG", 3)

# The product space
V = V0 * V1

# UFL mesh coordinates
x, y = SpatialCoordinate(mesh)

# Define and normalize cohomology generator
q = Function(V1).interpolate(
    as_vector([-y/(x**2+y**2), x/(x**2+y**2)]))
q /= sqrt(assemble(inner(q,q)*dx))

# Firedrake syntax for nullspace
V_basis = VectorSpaceBasis({q})
nullspace = MixedVectorSpaceBasis(V,[V.sub(0),V_basis])

# Function to decompose
f = as_vector([y+x, -x])

# Remove the contribution from cohomology
f_h = q * assemble(inner(f,q)*dx)
f_bar = f - f_h

# print("normalization check:")  # should be near 0
# print(assemble(inner(f_bar,q)*dx))

# Firedrake TnT functions
sigma, u = TrialFunctions(V)
tau, v = TestFunctions(V)

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

# The linear form
L = inner(f_bar, v) * dx

# These solver options are a complete guess
# Source https://github.com/ethrelfall/Finite-element-exterior-calculus/blob/main/source_problem_on_annulus.py
sp_it = {
   "ksp_type": "gmres",
   "pc_type": "fieldsplit",
   "pc_fieldsplit_type": "schur",
   "pc_fieldsplit_0_fields": "0",
   "pc_fieldsplit_1_fields": "1",
   "pc_fieldsplit_schur_precondition": "selfp",
   "fieldsplit_0_pc_type": "ilu",
   "fieldsplit_0_ksp_type": "preonly",
   "fieldsplit_1_ksp_type": "preonly",
   "fieldsplit_1_pc_type": "gamg",
   "fieldsplit_1_mg_levels_pc_type": "sor",
   "ksp_monitor": None,
#    "ksp_max_it": 500,
   "snes_monitor": None,
#    "ksp_rtol": 1.0e-10
}

# Solve the system
g = Function(V)
solve(a == L, g, nullspace=nullspace, solver_parameters=sp_it)

# Extract Hodge components and interpolate to V1 for plotting
sigma, u = g.subfunctions
grad_sigma = Function(V1).interpolate(grad(sigma))
delta_du = Function(V1).interpolate(f_bar-grad(sigma))
ff = Function(V1).interpolate(f)
f_harm = Function(V1).interpolate(f_h)

# Numerical checks
print(f'(dsigma, f_h): {assemble(inner(grad_sigma,f_h)*dx)}')
print(f'(deltadu, f_h): {assemble(inner(delta_du,f_h)*dx)}')
print(f'(deltadu, dsigma): {assemble(inner(delta_du,grad_sigma)*dx)}')
print(f'(u,q): {assemble(inner(u,q)*dx)}')

### Additional check that coexact part matches a direct solve

# Bilinear/Linear forms
w = TrialFunction(V1)
v = TestFunction(V1)
a = inner(w, v) * dx
# u-=q*(assemble(inner(u,q)*dx)) # Normalization
L = inner(curl(u), curl(v)) * dx

# Solve and print
d = Function(V1)
solve(a == L, d)
print(f'|delta_du-w|^2: {norm(d-delta_du)}')

# Plot the original vector field and its decomposition
fig, ax = plt.subplots(1, 4, figsize=(16,4))
# firedrake.tricontourf(sigma, axes=ax[0], cmap='cividis', levels=20)
firedrake.streamplot(ff, **{'resolution': 1/32}, axes=ax[0], cmap='cividis')
firedrake.streamplot(grad_sigma, **{'resolution': 1/32}, axes=ax[1], cmap='cividis')
firedrake.streamplot(delta_du, **{'resolution': 1/32}, axes=ax[2], cmap='cividis')
firedrake.streamplot(f_harm, **{'resolution': 1/32}, axes=ax[3], cmap='cividis')
# firedrake.quiver(f_harm, axes=ax[1], cmap='gray')
# firedrake.quiver(grad_sigma, axes=ax[0], cmap='gray')
# firedrake.quiver(delta_du, axes=ax[2], cmap='gray')
plt.show()

[dham]

My FEEC is certainly weaker than yours so I'm a bit out of my depth here. However I note that you have an iterative solver with the default tolerance. I can't remember what that is but that will introduce an error in g that presumably flow through into your hodge decomposition. Does decreasing ksp_rtol change the error in the orthogonality calculation?

[agrubertx]

Thanks for your response! The problem indeed appears to be something with the solver options that I do not understand. I have tried the following:

• Set ksp_rtol to 1.0e-11 (anything lower gave DIVERGED_LINEAR_SOLVE even though the printouts were converging...?) This did not work, and the orthogonality errors remained essentially unchanged.
• Change the FEM order of both spaces to 2 instead of 3, and set ksp_rtol to 1.0e-11. This worked, producing orthogonality errors on the order of $10^{-10},10^{-11}$, which is more like what I expect.
• With ksp_rtol at 1.0e-11 and FEM order at 2, perform a quick-and-dirty refinement study by cutting the mesh size in half 3 times. Everything converged as expected, with the weak $(u,q)=0$ condition becoming better and better satisfied.

Ultimately, I am more convinced now that my implementation is correct. I am still very interested in the two questions below the figure in my OP, though. Please let me know if you have thoughts about those as well.

Anthony

[dham]

I think in general doing 7.1 directly in finite element is hard because you'd need a finite element space that represents the space of harmonic forms, and that's not so straightforward. That's the reason that the approach of projecting the harmonic forms out of the solver is adopted.

For your second question, yes I think numerically constructing the VectorSpaceBasis is probably fine (especially because a lot of iterative solvers are reasonably null-space robust anyway, so approximately projecting out the null space is probably enough).

[agrubertx]

I see, thanks for your thoughts. I will close this thread now, but I may bug you again later with a new one as I try more things!

Best,
Anthony