Poisson equation the stubborn way

variational formulation
1

I am trying to understand what is a good variational formulation and what is a bad variational formulation and why. It seems common knowledge, that one should minimize the degree of the highest differential operator. Here is an example of a naive Poisson equation where I did not do that:

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
# F = inner(-grad(phi), -grad(v)) * dx - inner(rho, v) * dx # works
bc = DirichletBC(V, Constant(0.0), "on_boundary")
solve(F == 0, phi, bc)

It throws an error:

*** -------------------------------------------------------------------------
*** Error:   Unable to solve nonlinear system with NewtonSolver.
*** Reason:  Newton solver did not converge because maximum number of iterations reached.
*** Where:   This error was encountered inside NewtonSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------

It even throws an error, if I give it an analytic solution as initial guess:

phi = interpolate(Expression("-0.5*(x[0]-1)*(x[0]+1)", degree=2), V)
  • Why this error?
  • Why is inner(grad(phi), grad(v)) better then inner(-div(grad(phi)), v)?
  • Can I change the definition of the vector space V to make inner(-div(grad(phi)), v) work?
2

Someone else may answer your question specifically, but for a grasp of the fundamentals I recommend reading the introductory literature.

One of many which I like:
The Finite Element Method: Linear Static and Dynamic Finite Element Analysis By Thomas J. R. Hughes

Another I like:
Numerical Solution of Partial Differential Equations by the Finite Element Method by Claes Johnson

1 Like
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
4

Thanks this answer was quite enlightening. Is the fact that dx only integrates of the interiors of elements documented somewhere?

5

It comes down to terminology, since “integrates” is typically understood in the sense of Lebesgue when formulating finite element methods, so integrating over element interiors is equivalent to integrating over the full domain, and it’s simply invalid to talk about integrating an object like -\nabla\cdot(\nabla\phi^h), since it’s not a function. (It’s arguably an abuse of notation to write things like \int \delta_0(x)v(x)\,dx = v(0); really it’s the application of a linear functional \delta_0 to the function v, i.e., \langle \delta_0, v\rangle.) From that perspective, it’s a user error to ask for an ill-defined integral. However, that’s not the full story, since many residual-based stabilized methods actually do include terms that are only integrable on element interiors and should technically be written using a sum over element integrals, which is one reason I can think of why things like inner(-div(grad(phi)),v)*dx are allowed to compile without errors.

EDIT: To answer the question more directly, yes, the UFL manual explicitly defines the dx measure as a sum over element integrals. See the formula at the bottom of page 4.

1 Like