Maxwell Discretizations: The Good, The Bad & The Ugly¶
Jay Gopalakrishnan
MTH 653 Class Notes, Spring 2019
Are Nedelec elements really necessary to solve Maxwell equations? Imagine a world without the Nedelec element, however terrible that might be: a world where people were forced to solve Maxwell equations using the Lagrange element. How would one go about using Lagrange elements? Does one get a reasonable method? We shall answer all these questions in this notebook. Here are the punchlines in advance.
The Good¶
Compare the $x$-component of an exact electric field (left) with its approximation computed using the Nedelec element (right):
The Bad¶
Compare the $x$-component of an exact electric field (left) with its approximation computed using the Lagrange element (right):
The Ugly¶
As mesh is made finer the Lagrange discretization converges at a reasonable rate, but it does not converge to the exact solution! In numerical analysis, this is one of the worst kind of failures that one can expect.
We now proceed to see all this in detail.
Geometry¶
The problems we intend to exemplify are only visible in a non-convex domain, so we construct an L-shaped domain.
import ngsolve as ng
import netgen.gui
%gui tk
from netgen.geom2d import SplineGeometry
from ngsolve import CoefficientFunction, InnerProduct, Integrate
from ngsolve import BilinearForm, LinearForm, dx, Mesh, Draw, HCurl, H1
from ngsolve import GridFunction
from ngsolve import sin, cos, asin, acos, x, y, z, curl, grad, atan2, sqrt
from prettytable import PrettyTable
from math import pi
import numpy as np
#
# (-1,1) G6 (1,1)
# +----------------------+
# | |
# | | G4
# | (0,0) |
# G3 | +----------+ (1,0)
# | | G2
# | |G1
# | G5 |
# +-----------+
# (-1,-1) (0,-1)
#
geo = SplineGeometry()
ptlist = [(0,0), (1,0), (1,1), (-1,1), (-1,-1), (0,-1)]
pts = [geo.AppendPoint(*p) for p in ptlist]
geo.Append(['line', pts[0], pts[1]], bc='midhoriz', leftdomain=1, rightdomain=0)
geo.Append(['line', pts[1], pts[2]], bc='rghtvert', leftdomain=1, rightdomain=0)
geo.Append(['line', pts[2], pts[3]], bc='tophoriz', leftdomain=1, rightdomain=0)
geo.Append(['line', pts[3], pts[4]], bc='leftvert', leftdomain=1, rightdomain=0)
geo.Append(['line', pts[4], pts[5]], bc='bothoriz', leftdomain=1, rightdomain=0)
geo.Append(['line', pts[5], pts[0]], bc='midvert', leftdomain=1, rightdomain=0)
Boundary value problem¶
Let $\Omega$ be the L-shaped domain we constructed above. The problem is to find $E$ satisfying $$\DeclareMathOperator{\curl}{curl}\DeclareMathOperator{\div}{div} \DeclareMathOperator{\grad}{grad} \begin{aligned} \curl \curl E & = 0 && \text{ in } \Omega \\ \div E & = 0 && \text{ in } \Omega \\ E\cdot t & = g && \text{ on } \partial\Omega. \end{aligned} $$ The data $g$ is assumed to have an extension $E_g \in H(\curl) \cap H(\div)$ such that $ E_g \cdot t = g $ on $\partial\Omega$, so that problem of finding $E$ reduces to the problem of finding $E_0 = E - E_g$ with homogeneous essential boundary conditions.
Make do with Lagrange elements¶
One might be tempted to argue that controlling $\curl$ and $\div$ essentially controls all first order derivatives, so the Maxwell solution might satisfy an $H^1$-formulation. Indeed, writing $E = E_0 + E_g$, the we immediately see that the unknown function $E_0$ belongs to the space $$ X = H_0(\curl) \cap H(\div). $$ This results in the following weak formulation: Find $E = E_0 + E_g$ such that $E_0 \in X$ satisfies $$ \begin{aligned} (\curl E, \curl v) + (\div E, \div v) & = 0, && \text{ for all } v \in X. \end{aligned} $$ A piecewise polynomial function is in $H(\curl)$ if its tangential component is continuous. It is in $H(\div)$ if its normal component is continuous across element interfaces. Hence, it is in $X$ if all its components are continuous across element interfaces. Thus, the Lagrange finite element space is a conforming finite element subspace of $X$. This formulation therefore appears to offer an avenue to compute the Maxwell solution using the commonly available Lagrange finite elements.
The standard conforming method¶
One can also, as we have already seen previously, construct a weak formulation treating the divergence equation weakly: Find $E$ in $H(\curl)$, $E = E_0 + E_g,$ where $E_0 \in H_0(\curl)$ satisfies, together with $\phi \in H_0^1$, $$ \begin{aligned} (\curl E, \curl v) + (\grad \phi, v) & = 0 && \text{ for all } v \in H_0(\curl) \\ (\grad\psi, E) & = 0 && \text{ for all } \psi \in H_0^1. \end{aligned} $$ For this formulation, we use the Nedelec space for approximating $E$ and the Lagrange finite element space for approximating $\phi$.
We shall see below that widely different approximations are obtained using the above two approaches.
Choosing an exact solution¶
To see the above-mentioned effects, we must choose an exact solution of low regularity. The exact solution we shall use, in polar coordinates, is $$ E= \grad\left(r^{2/3} \sin(2\theta/3)\right). $$ Note the following regarding $E$:
- $E$ satisfies the above boundary value problem with a $g$ that is nonzero only on the edges of $\Omega$ that do not meet the origin.
- Although $E$ has a singularity, the data $g$ is smooth along $\partial\Omega$, hence it offers a case where no special integration is required to assemble the corect right hand side.
- The components of $E$ are not in $H^1(\Omega)$.
We need to implement this function for error computation. For this, we need $\theta$, which can be obtained using arcsin, arccos or arctan. But please do be careful: for example, see the result of setting theta = acos(x/r):
r = sqrt(x*x + y*y)
theta = acos(x/r)
mesh = Mesh(geo.GenerateMesh(maxh=1/8))
Draw(theta, mesh, 'theta')
The values are incorrect
- in the third quadrant,
- and near the origin due to division by (close to) zero.
These problems can be solved as follows:
r = sqrt(x*x + y*y)
rinv = 1.0/ng.IfPos(r-1e-15, r, 1.e-15) # threshold to avoid 0-division
theta = ng.IfPos(y+1e-15, ng.acos(x * rinv), pi - ng.asin(y*rinv))
Draw(theta, mesh, 'theta')
Alternately, one may use arctangent with two arguments (y and x, avoiding division by zero).
theta = ng.atan2(y, x)
Draw(theta, mesh, 'theta')
But again here, we must be careful, as it produces a branch cut through our $\Omega$. We rotate the coordinate system, take arctan, and rotate back to put the branch cut outside our domain.
alpha = -(pi/2 + pi/4)
rotatedx = sin(alpha) * x - cos(alpha) * y
rotatedy = sin(alpha) * x + cos(alpha) * y
theta = ng.atan2(rotatedy, rotatedx) - alpha
Draw(theta, mesh, 'theta')
Using either of these "fixed up" $\theta$, we proceed to define the exact solution:
Eexact = ((2/3)*(cos(theta)*sin(2*theta/3) -
sin(theta)*cos(2*theta/3))*pow(rinv,1/3),
(2/3)*(cos(theta)*cos(2*theta/3) +
sin(theta)*sin(2*theta/3))*pow(rinv,1/3))
Eexact = ng.CoefficientFunction(Eexact)
Draw(Eexact, mesh, 'Eexact')
The Nedelec approximation¶
def SolveByNedelec(mesh, bc_E, bc_phi, p=1):
""" Given boundary data for E.t in bc_E and for phi in bc_phi, solve
by the above-mentioned Nedelec approach.
"""
V = HCurl(mesh, type1=True, order=p, dirichlet='[a-z]*')
L = H1(mesh, order=p, dirichlet='[a-z]*')
X = ng.FESpace([V, L])
u, phi = X.TrialFunction()
v, psi = X.TestFunction()
a = BilinearForm(X, symmetric=True)
a += (curl(u)*curl(v) + grad(phi)*v + u*grad(psi)) * dx
f = LinearForm(X)
Ephi = GridFunction(X, 'E_Nedelec')
with ng.TaskManager():
Ephi.components[0].Set(bc_E, ng.BND)
Ephi.components[1].Set(bc_phi, ng.BND)
a.Assemble()
f.Assemble()
r = f.vec.CreateVector()
r.data = f.vec - a.mat * Ephi.vec
Ephi.vec.data += a.mat.Inverse(X.FreeDofs()) * r
return Ephi, X
We apply the above routine by giving the boundary condition argument as the exact solution. While Eexact is the exact E-component of the solution, the exact phi is $0$.
Ephi, X = SolveByNedelec(mesh, Eexact, CoefficientFunction(0), p=4)
Draw(Ephi.components[0], mesh, 'ENedelec')
The solution picture should look very similar to what was described in "The Good" punchline above.
The Lagrange approximation¶
def SolveByLagrange(mesh, bc_E, p=1):
""" Solve using the above-mentioned approach using Lagrange
elements (only) for each component of the electric field,
given E.t boundary data in bc_E. """
# Make Lagrange spaces so that their product has the
# required tangential boundary conditions:
Vx = H1(mesh, order=p, dirichlet='midhoriz|tophoriz|bothoriz')
Vy = H1(mesh, order=p, dirichlet='rghtvert|leftvert|midvert')
X = ng.FESpace([Vx, Vy])
ux, uy = X.TrialFunction()
vx, vy = X.TestFunction()
u = CoefficientFunction((ux, uy))
v = CoefficientFunction((vx, vy))
# Define the two differential operations required for the form:
def curl2D(w0, w1):
dw0 = grad(w0)
dw1 = grad(w1)
return dw1[0] - dw0[1]
def div2D(w0, w1):
dw0 = grad(w0)
dw1 = grad(w1)
return dw1[1] + dw0[0]
# System:
a = BilinearForm(X, symmetric=True)
a += (curl2D(ux, uy) * curl2D(vx, vy) + \
div2D(ux, uy) * div2D(vx, vy)) * dx
f = LinearForm(X)
# Solve:
u = GridFunction(X, 'E_Lagrange')
u.components[0].Set(bc_E[0], ng.BND)
u.components[1].Set(bc_E[1], ng.BND)
a.Assemble()
f.Assemble()
r = f.vec.CreateVector()
r.data = f.vec - a.mat * u.vec
u.vec.data += a.mat.Inverse(X.FreeDofs()) * r
return u, X
mesh = ng.Mesh(geo.GenerateMesh(maxh=1/8))
E, X = SolveByLagrange(mesh, Eexact, p=4)
Draw(E.components[0], mesh, 'ELagrange')
This is clearly the solution shown in "The Bad" punchline above. Its values near the non-convex corner shows smooth variations, indicating that it has completely missed the singularity!
Convergence study¶
You might wonder if this situation gets remedied on finer meshes. To study this, we perform a convergence study: Start with a (coarse) mesh, and solve the same Maxwell problem on successively refined meshes. Each refinement below is obtained by connecting the midpoints of edges of each triangle in the current mesh.
Rates for Lagrange approximation¶
def SolveByLagrangeSuccessive(hcoarse=1/4, p=1, nrefinements=5):
"""
Starting with a mesh of grid size "hcoarse", solve using the
Lagrange method on successively refined meshes. Store all solutions,
spaces and meshes and return them in lists.
"""
Es = []
Xs = []
meshes = []
mesh = ng.Mesh(geo.GenerateMesh(maxh=hcoarse))
for ref in range(nrefinements):
meshes.append(ng.Mesh(mesh.ngmesh.Copy()))
u, X = SolveByLagrange(mesh, Eexact, p=p)
Es.append(u)
Xs.append(X)
mesh = meshes[-1]
mesh.ngmesh.Refine()
return Es, Xs, meshes
Given solutions $E_i$ on a sequence of meshes of grid size $h_i= h_0/2^i$, we can estimate the rate of convergence by examining at what rate $$ \| E_i - E_{fine} \|_{L^2(\Omega)} \to 0. $$ Here $E_{fine}$ is the numerical solution computed on the finest mesh (max $i$). This is usually how we compute the NOC (Numerical Order of Convergence) of a method when we have no access to the exact solution.
In this problem however, we know the exact solution, so we are in a position to cross check the NOC with how the exact error $$ \| E_i - E_{exact} \|_{L^2(\Omega)} \to 0. $$
def convergence_study(Es, Xs, meshes):
"""Given solutions on successively refined meshes, return
|| E_i - E_fine|| and ||E_i - E_exact|| for i-th mesh, for all i.
"""
E_diff = []
E_err = []
fine_E = ng.GridFunction(Xs[-1])
with ng.TaskManager():
fine_E.components[0].Set(Es[-1].components[0])
fine_E.components[1].Set(Es[-1].components[1])
for i in range(len(meshes)-1):
diffE0 = fine_E.components[0] - Es[i].components[0]
diffE1 = fine_E.components[1] - Es[i].components[1]
dE = CoefficientFunction((diffE0, diffE1))
E_diff.append(np.sqrt(Integrate(dE*dE, meshes[i])))
dE = CoefficientFunction((Eexact[0] - Es[i].components[0],
Eexact[1] - Es[i].components[1]))
E_err.append(np.sqrt(Integrate(dE*dE, meshes[i])))
return np.array(E_diff), np.array(E_err)
Es, Xs, meshes = SolveByLagrangeSuccessive(hcoarse=1/4, p=1,
nrefinements=6)
E_diff, E_err = convergence_study(Es, Xs, meshes)
E_diff
These numbers certainly look like they are converging. Let's make a quick function to print the rate of convergence:
def tabrate(name, dat):
col = ['h', name, 'rate']
t = PrettyTable()
t.add_column(col[0], ['1/'+str(2**(2+i)) for i in range(len(dat))])
t.add_column(col[1], ['%.7f'%e for e in dat])
t.add_column(col[2], ['*'] + \
['%1.2f'%r for r in np.log(dat[:-1]/dat[1:])/np.log(2)])
print(t)
tabrate('||Eh-Efine||', E_diff)
This seems to indicate that the method converges.
The problem is only revealed when we see the exact errors:
tabrate('||Eh-Eexact||', E_err)
E_err
These errors do not go to zero, thus showing that the method converges, but to something distant from the exact solution!
Rates for Nedelec approximation¶
We now repeat the same study as above, but now using the "The Good" method with Nedelec elements.
def SolveByNedelecSuccessive(hcoarse=1/4, p=1, nrefinements=6):
Ephis = []
Xs = []
meshes = []
mesh = ng.Mesh(geo.GenerateMesh(maxh=hcoarse))
for ref in range(nrefinements):
meshes.append(ng.Mesh(mesh.ngmesh.Copy()))
Ephi, X = SolveByNedelec(mesh, Eexact, CoefficientFunction(0), p=p)
Ephis.append(Ephi)
Xs.append(X)
mesh = meshes[-1]
mesh.ngmesh.Refine()
return Ephis, Xs, meshes
def convergence_study2(Ephis, Xs, meshes):
E_diff = []
E_err = []
fine_Ephi = ng.GridFunction(Xs[-1])
Efine = fine_Ephi.components[0]
with ng.TaskManager():
Efine.Set(Ephis[-1].components[0])
for i in range(len(meshes)-1):
dE = Efine - Ephis[i].components[0]
E_diff.append(sqrt(Integrate(dE*dE, meshes[i])))
dE = Eexact - Ephis[i].components[0]
E_err.append(sqrt(Integrate(dE*dE, meshes[i])))
return np.array(E_diff), np.array(E_err)
Ephis, Xs, meshes = SolveByNedelecSuccessive(hcoarse=1/4, p=1,
nrefinements=6)
Ediff, Eerr = convergence_study2(Ephis, Xs, meshes)
tabrate('||Eh-Efine||', Ediff)
tabrate('||Eh-Eexact||', Eerr)
Clearly the exact errors and the numerically estimated errors are converging to zero.
In the accompanying lecture, we shall the mathematics behind the "The Bad" alternative and why it fails. To my knowledge, the possibility of this failure was first hinted at in a theoretical paper by Martin Costabel. We have already shown the nice convergence properties of "The Good" alternative in prior lectures. (Can you explain the rate of 0.66 that we see above from the theory in class?)