Poisson equation the stubborn way

variational formulation
3

I second Nate’s suggestion to take some time with the fundamentals. To answer the questions specifically,

  • Why this error? Because dx is really a sum of integrals over interiors of elements. You can see that, on the interior of each element, -\nabla\cdot(\nabla\phi) is zero when \phi is piecewise linear. Thus, the left-hand side matrix is singular.

  • Why is inner(grad(phi), grad(v)) better? Because the other definition is not correct! Glossing over some technicalities about the theory of distributions, recall that the derivative of a Heaviside function is a Dirac delta. Considering -div(grad(phi)) on the element interiors only (as per the dx measure) is missing the Dirac deltas that should be part of the derivative of the discontinuous function grad(phi). See the following example, where dS is a sum of integrals over interior mesh facets:

from fenics import *

mesh = IntervalMesh(10, -1.0, 1.0)
V = FunctionSpace(mesh, "P", 1)
phi = Function(V)
rho = Constant(1.0)
v = TestFunction(V)
F = inner(-div(grad(phi)), v) * dx - inner(rho, v) * dx # does not work

# Because grad(phi) is discontinuous, its derivative includes Dirac deltas
# which must be explicitly added, because *dx is a sum of integrals
# over element interiors.
n = FacetNormal(mesh)
F += inner(jump(grad(phi),n),avg(v))*dS

bc = DirichletBC(V, Constant(0.0), "on_boundary")
solve(F == 0, phi, bc)

# Compare with standard formulation:
phi2 = Function(V)
F = inner(-grad(phi2), -grad(v)) * dx - inner(rho, v) * dx # works
bc = DirichletBC(V, Constant(0.0), "on_boundary")
solve(F == 0, phi2, bc)

# Solutions are identical:
phiDiff = phi - phi2
print(assemble((phiDiff**2)*dx))
  • Can I change the definition of the vector space V to make inner(-div(grad(phi)), v) work? Mathematically, yes, you could use a C^1 space of functions, in which case the jump terms in the above example vanish, although this is not supported directly in FEniCS. (You can do it using a library I worked on, though.) However, it would then be exactly equivalent to inner(grad(phi),grad(v)) (with Dirichlet BCs).
5 Likes