Time Harmonic Maxwell equation - meshing influences solutions

variational formulation
1

I am trying to solve the vectorial helmholtz equation in a cylindrical inhomogenous waveguide, with a geometry and olverlayn permittivity shown in Fig. (a). The variational formulation for the problem with PEC boundary conditions and E split into a transverse and z term is shown in (b) and its implementation for Fenics is based on the work from Lezar. A eigenvalue spectrum is shown in (c), with the horiztonal line showing the maximum physically possible. I guess the higher values belong to spurious modes, fine. (d) shows the x-component of a mode. For some mesh parameters I get solutions that make sense (as given in the full code below). However, when making the mesh finer by adjusting

mesh = generate_mesh(domain, 15)

more spurious modes kick in and my previous solutions are overlayn with random new features (e). I haven’t seen this kind of behaviour in COMSOL before. Have I formulated an ill-posed problem maybe? Using gmsh to refine the mesh only around the core also did not help. A ‘minimal’ example is given below.

This has been bugging me for a quite a while, so if someone knows a way out I would appreciate it.
Steffen

text4535-1
text4535-12456×2152 534 KB

from dolfin import * 
import numpy as np
import fenics as fn
from mshr import *

class set_properties(UserExpression):
    def __init__(self, subdomains, param, **kwargs):
        super().__init__(**kwargs)
        self.subdomains = subdomains
        self.params = param

    def eval_cell(self, values, x, cell):
        values[0] = self.params[self.subdomains[cell.index]-1]
        
class OperatingWavenumberSquared(UserExpression):
    def eval (self, values, x):
        values[0] = self.__k_o_squared
    def set_wavelength(self, wavelength):
        self.__k_o_squared = (2*np.pi/wavelength)**2
        
def curl_t(w):
    return Dx(w[1], 0) - Dx(w[0], 1)

vector_order = 2
nodal_order = 3

wavelength = 1.55
k0 = 2*pi/wavelength

r_clad = 40
r_core = 8

n_clad = 1.4378
n_core = 1.4457

domain = Circle(Point(0, 0), r_clad)
core = Circle(Point(0, 0), r_core)

domain.set_subdomain(1, core)
mesh = generate_mesh(domain, 15)

surface_markers = MeshFunction('size_t', mesh, 2, mesh.domains())
dx = Measure('dx', domain=mesh, subdomain_data=surface_markers)

e_rs = np.array([n_core**2, n_clad**1])
e_r = set_properties(surface_markers, e_rs, degree=2)

k_0_squared = OperatingWavenumberSquared()
k_0_squared.set_wavelength(wavelength)
> 
vector_element = FiniteElement("N1curl", mesh.ufl_cell(), vector_order)
nodal_element = FiniteElement("Lagrange", mesh.ufl_cell(), nodal_order)

V = FunctionSpace(mesh, vector_element * nodal_element)   # NOTE: V.sub(0) and V.sub(1) are the subspaces 

(N_i, L_i) = TestFunctions(V)
(N_j, L_j) = TrialFunctions(V)

u_r = 1.0
s_tt_ij = 1.0/u_r * dot(curl_t(N_i), curl_t(N_j))
t_tt_ij = e_r * dot(N_i, N_j)
s_zz_ij = 1.0/u_r * dot(grad(L_i), grad(L_j))
t_zz_ij = e_r * L_i * L_j

b_tt_ij = 1.0/u_r * dot(N_i, N_j)
b_tz_ij = 1.0/u_r * dot(N_i, grad(L_j))
b_zt_ij = 1.0/u_r * dot(grad(L_i), N_j)

a_ij = ( s_tt_ij - k_0_squared * t_tt_ij ) * dx
b_ij = ( s_zz_ij - k_0_squared * t_zz_ij + b_tt_ij + b_tz_ij + b_zt_ij ) * dx

electric_wall = DirichletBC(V, Constant((0.0, 0.0, 0.0)), DomainBoundary())

parameters["linear_algebra_backend"] = "PETSc"

A = PETScMatrix()
B = PETScMatrix()

assemble(a_ij, tensor=A)
assemble(b_ij, tensor=B)

electric_wall.apply(B)
electric_wall.apply(A)

neigs = 50
disp = np.zeros(neigs)

# Create VTK file for saving solution
vtkfile = File('solved_modes/solution.pvd')
    
solver = SLEPcEigenSolver(A, B)

solver.parameters['tolerance']          = 1e-12
solver.parameters["spectrum"]           = "target magnitude"
solver.parameters["spectral_transform"] = "shift-and-invert"
solver.parameters["spectral_shift"]     = -n_core**2*k0**2
solver.parameters["solver"]             = "arnoldi"

solver.solve(neigs)

assert solver.get_number_converged() > 0

for j in range(min(neigs, solver.get_number_converged())):
    r, c, rx, cx = solver.get_eigenpair(j)

    disp[j] = np.sqrt(np.abs(r/k0**2))

    mode = Function(V)

    mode.vector()[:] = rx

    (Et, Ez) = mode.split()
    Et.rename("Et", "Transverse mode solution")
    Ez.rename("Ez", "z-axis mode solution")
    vtkfile << (Et, j)
    vtkfile << (Ez, j)
2

You write “Helmholtz”, but do you mean “Maxwell”?

Are you sure this isn’t a visualisation issue? Writing the vertex values of a Nedelec function for visualisation is iladvised. Perhaps project to a C^1 function space since your polynomial approximation of the electric field is degree 2.

Regarding the spurious modes: Won’t these just arise due to the boundary conditions you’ve implemented? If you’re using H(\mathrm{curl})-conforming elements you shouldn’t see other spurious modes. Maybe I’m wrong here?

For example: you could remove the rows and columns corresponding to the DirichletBCs from A and B using petsc4py. Or set the diagonal values corresponding to the DirichletBC DoFs to some large number (DOLFIN’s default is 1.0) so not to affect the modes close to your spectral shift.

One more sanity check: Do the trial and test functions L_i and L_j correspond to your construction of an involution term in the richer space? May be worth checking the corresponding FE solution is 0 everywhere (if I remember the analysis correctly).

3

thanks for your response.

maybe vector Helmholtz equation is not a thing. I should have better said, I am trying to solve the time harmonic wave-equation (derived from Maxwell’s equation without sources)

\nabla \times \mu^{-1} \nabla \times \vec{E} - \omega^2 \epsilon \vec{E} = 0

with PEC boundary conditions.

What I was showing in d, e was vector solution along x. How would I do the projection you are proposing in Fenics?

(Et, Ez) = mode.split()
xhat = fn.Expression((‘1.0’, ‘0.0’), V.sub(0), degree=1)
Ex = fn.dot(xhat, Et)

I also thought that was what the shift-and-invert was for to shift your solutions away from the spurious modes. But then I saw these eigenvalues that shouldn’t be there and I cannot explain to myself why they are there.

Can you explain what you mean with:

I will check your other suggestion about the diagonal values of the DirichletBC DoFs.

4

I’m pretty rusty in this area. Hopefully an expert appears soon. Nevertheless:

I was wondering about the involution term which is typically used to enforce Gauß’s law. Perhaps this isn’t necessary in the generalised eigenvalue problem.

To do the projection you can just use FEniCS’s project function. E.g. E_cg = project(E, VectorFunctionSpace(mesh, "CG", 1)). You’ll get the best approximation of E in a CG1 space according to CĂ©a’s Lemma. This should then give you a better function for visualisation with the output method you’re using to VTKFile.

Shift-and-invert alters the nature of the eigenvalue problem so that eigenvalues close to the spectral shift are found first by the iterative solution method due to inversion about that point making them the largest in a modified problem. This does not prevent spurious modes from appearing. My understanding of a spurious mode is that it constitutes a numerical approximation which insidiously appears to converge to a solution. But that solution is not a true solution of the original Maxwell problem. See for example this excellent demo.

The modes that appear due to the DirichletBCs correspond to the zeroed rows and columns with a single diagonal entry corresponding to that degree of freedom. The simple workaround for this is to set the diagonal entries to be much larger than the unity default set by dolfin. This means that when you perform the spectral shift, the eigenvalue solutions corresponding to the DirichletBCs will be very far away from your shift target.

2 Likes
5

I tried two more things
First:the projection, while the output looks a little better (right) the real problem seems to be somewhere else.

I used

r, c, rx, cx = solver.get_eigenpair(5)
mode = Function(V)
mode.vector()[:] = rx
(Et, Ez) = mode.split()
E_cg = project(Et, VectorFunctionSpace(mesh, “CG”, 1))
plot(E_cg)

projection_to_CG_space
projection_to_CG_space1913×900 228 KB

Second: Setting the diagonal values of my matrices corresponding to the DirichletBCs gives me trouble.
I successfully located the positions that I should change in my matrices through

electric_wall = DirichletBC(V, Constant((0.0, 0.0, 0.0)), DomainBoundary())
boundary_points = np.asarray(list(electric_wall.get_boundary_values().keys()))

but an attempt to change the values directly in the matrices through A[mask] = some_value didn’t cut it. Is there a different approach how to change the values that DirichletBC places when applying the bc’s?

6

It has been a while since I’ve looked at this particular problem, but I recall that the mode spectrum for lossless inhomogeneous waveguides is complex. You will need to change your formulation to take this into account, as complex eigenmodes appear with eigenvalues in complex pairs (the eigenvalue \gamma^2 will be complex). Formulation errors are likely polluting your solution.

Take a look at
K. Zaki and Chunming Chen, “Complex modes in dielectric loaded waveguides,” 1987 Antennas and Propagation Society International Symposium, Blacksburg, VA, USA, 1987, pp. 8-11.

7

@BillS Thank you! That sounds like a good suggestion. I followed a previous comment from @nate and started reformulating my problem. I got stuck with the sesquilinearity of my variational formulation though, (u, v) = \int_\Omega u \overline{v} \; \mathrm{d} x = (u_r + j u_j, v_r - j v_j). Do I use this to flip signs when I define the elements of my final matrices (in particular the ‘imaginary’ system)? For example here,

s_tt_i = 1.0/u_r * inner(curl_t(v_N_i), curl_t(u_N_i))
t_tt_i = e_r * inner(v_N_i, u_N_i)
s_zz_i = 1.0/u_r * inner(grad(v_L_i), grad(u_L_i))

Creating a new function for inner wouldn’t really cut it no?

The other changes I implemented were:

  • mixed function space with both real and imaginary part (where N denote vector and L Lagrange elements)

    W = FunctionSpace(mesh, MixedElement([vector_element, nodal_element, vector_element, nodal_element]))
    U, V = TrialFunction(W), TestFunction(W)

    u_N_r, u_L_r, u_N_i, u_L_i = split(U)
    v_N_r, v_L_r, v_N_i, v_L_i = split(V)

  • write out my original formulation for both real and imaginary, assemble matrices, apply bcs add both together and solve.

    electric_wall = DirichletBC(W, Constant((0.0, 0.0, 0.0, 0.0, 0.0, 0.0)), DomainBoundary())
    A_i = PETScMatrix()
    

    assemble(a_i, tensor=A_i)
    

    electric_wall.apply(A_i)
    

    A = A_r + A_i

8

Yes, your matrix will look like

\int\limits_\Omega u\overline{v} dx = \langle v_r, u_r\rangle + \langle v_j, u_j\rangle + j\left (\langle v_r, u_j\rangle - \langle v_j, u_r\rangle\right )

the way you have it written. The resulting matrix will always have this type of block skew-symmetric form.

9

@BillS Great, what I don’t understand is if I start with such initial problem

\nabla \times \mu^{-1} \nabla \times \vec{E} - \omega^2 \epsilon \vec{E} = 0
\vec{n} \times \vec{E} = f~\mbox{ on $\Gamma_E$}

and say that E has both imaginary and real part E = E^r + i E^i . This however, leads to my initial system

\nabla \times \mu^{-1} \nabla \times \vec{E^r} - \omega^2 \epsilon \vec{E^r} = 0
\nabla \times \mu^{-1} \nabla \times \vec{E^i} - \omega^2 \epsilon \vec{E^i} = 0
\vec{n} \times \vec{E^r} = \vec{n} \times \vec{E^i} = f~\mbox{ on $\Gamma_E$}

So I would just end up with my initial variational formulation, this time for both real and imaginary part. Looking at the inner products both arguments are always real, so I wouldn’t even need a sesquilinear inner product. What am I missing here?

1 Like
10

You are absolutely correct. This is because the real and imaginary parts of the field must satisfy the same wave equation. What are the ways to get coupling between real and imaginary parts of the field?

  1. If there are lossy materials present. This is easy to see. But since we are dealing with the lossless problem, we need to look for other means.
  2. Through boundary conditions. A port or radiation boundary condition looks like a resistive load, hence will couple the real and imaginary field components. (I have a post in this in this discussion showing this boundary condition). This is the case for scattering that you show in your last post
  3. The other way is in situations where hybrid modes are present (not pure TE or TM). This is your waveguide case. This happens because the boundary condition between materials sometimes requires the existence of complex modes.

Your waveguide eigenvalue system is indefinite, which means that under certain circumstances it may yield complex conjugate pair eigenvalues. This will happen in certain frequency ranges for certain geometries. Also, \gamma^2 is complex. Your assembled matrices and DoF vectors must reflect this fact when you construct them.

2 Likes
11

@BillS Thank you for the clarification! I see that my boundary conditions and real materials do not lead to coupling between real and imaginary parts of the field. I checked my notes again and yes you are right, the modes that we expect here are hybrid indeed. But if I understand correctly this does not influence my choice of finite elements itself (N x CG).

For the formulation I got stuck. Starting from my original problem

\begin{array}{l} \nabla \times\left( \mu_r^{-1} \nabla \times E\right) - k_0^2 \epsilon_r E = 0 \end{array}

and having splitting E into transverse and axial part (\vec{E_t} + E_z \hat{z}) exp(-\gamma j z) I end up here

\begin{array}{l} \nabla_t \times \left( \mu_r^{-1} \times E_t\right) - \mu_r^{-1} \left( j \gamma\nabla_t E_z - \gamma^2 E_t\right) = k_0^2 \epsilon_r E_t \\ - \mu_r^{-1} \left( \nabla_t \cdot \left(\nabla_t E_z + j \gamma E_t\right) \right) = k_0^2 \epsilon_r E_z \end{array}

Are you suggesting I should split up \beta = \beta_r + j \beta_i and E = E_r + j E_i. Doing this I don’t see however how I would get towards a final form of

S \begin{bmatrix} E_{tr} \\ E_{zr} \\ E_{ti} \\ E_{zi} \end{bmatrix} = -\gamma^2 T \begin{bmatrix} E_{tr} \\ E_{zr} \\ E_{ti} \\ E_{zi} \end{bmatrix}

to be able to solve the system together?

12

I have not tried it myself, but on looking at the Dolfin eigenvalue demo, it looks like you can set up the system as
S\left [\begin{array}{c}E_t\\ E_z\end{array}\right ]-\gamma^2 T\left [\begin{array}{c}E_t\\ E_z\end{array}\right ]=0,
where \gamma^2, E_t, E_z are assumed complex. The matrices S and T are real. As far as I understand, you do not need to explicitly split the solution vectors. The SLEPcEigenSolver eigensolver returns complex data types.

Looking at your first post, the eigensolution looks correctly implemented. One thing I noted though: if you used the formulation shown in your first post, your E_t is really \gamma E_t. Did you take this into account?
In case you haven’t seen it, here is a link to the pages in Pelosi’s book describing the derivation.

13

@BillS ha, finally understood it. As my matrices S and T are not complex I don’t really have to split up my system into real and imaginary part. SLEPc can still give me complex solutions, so all should be fine. I erroneously fell into the ‘FENics cannot handle complex systems’ and had tunnel vision, in any case I should check out dolfinx at some point. I don’t have access to my sim-computer right now, but I will check the old code with sesquilinear inner product and that should hopefully do it.

14

@BillS @nate: I really though that sesquilinearity was my issue, but I saw that the form language has a sesquilinear inner product by default: Form language docs. Or is there a switch in SLEPc to enable complex results? Up to now my spectrum is completely real, with values obtained through

r, c, rx, cx = solver.get_eigenpair(j)

One thing I noted though: if you used the formulation shown in your first post, your Et is really ÎłEt . Did you > take this into account?
In case you haven’t seen it, here is a link 1 to the pages in Pelosi’s book describing the derivation.

Thank you for the book excerpt. You are right in the final solution I should scale my resulting Et vectors, haven’t done that yet.

I will try again to scale the DirichletBCs diagonal elements as noted by nate

The modes that appear due to the DirichletBCs correspond to the zeroed rows and columns with a single diagonal entry corresponding to that degree of freedom. The simple workaround for this is to set the diagonal entries to be much larger than the unity default set by dolfin. This means that when you perform the spectral shift, the eigenvalue solutions corresponding to the DirichletBCs will be very far away from your shift target.

15

Have you thought about taking things much more slowly and ensuring you’re solving the problem you want to solve? The first problems I’d solve to make sure what I’m doing is correct follow:

  1. Linear problem, find E_h \in V^h such that:

    (\nabla \times E_h, \nabla \times F) - \gamma^2 (E_h, F) = (f, F) \quad \forall F \in V^h

    where we prescribe the true solution E = (\sin(\gamma x), \sin(\gamma y))^\top (set BCs accordingly), f=0 and \gamma^2 is not an eigenvalue of the linear system. You can then compare the approximation E_h with E and make sure you’re converging at optimal rates.

  2. Eigenpair problem, find (E_{h,mn}, \lambda_{h,mn}\neq0) \in V^h \times \mathbb{R} such that:

    (\nabla \times E_{h,mn}, \nabla \times F) = \lambda_{h,mn} (E_{h,mn}, F) \quad \forall F \in V^h

    subject to appropriate BCs (n \times E_{h,mn} = 0 on \partial\Omega). I’m rusty here but I vaguely recall the true solution \nabla \times E_{mn} = \sin(m \pi x) \sin (n \pi y) and \lambda^2_{mn} = (m \pi)^2 + (n \pi)^2 when \Omega = (0, 1)^2. A quick google search for “Rectangular TM waveguide modes” should help.

Both these problems are real valued. Once I was happy with these I’d move on to a complex or more complicated formulation.

I’m writing these from memory so they could be filled with mistakes, but hopefully useful as a guide. I believe I recall these examples from Peter Monk’s excellent book.

16

Thank you for the great advice. Initially I did some testing with the eigenpair problem you mentioned in 2 and my obtained eigenvalues matched the theoretical values \lambda_{mn}^2 = (m\pi)^2 + (n \pi)^2. But I never looked into the eigenmodes and convergence too much, I will catch up on that
and have a deep look into your cited reference.

17

@nate I’ve analyzed the eigenpair problem, considering only the transverse field. For given problem with Dirichlet BCs

\nabla_t \times \left( \frac{1}{\mu_r} \nabla_t \times u_t\right) - k_c^2 \epsilon_r u_t = 0

It seems to work well, is the optimal convergence something you know from experience or is there something you can deduce from the element type used?

Capture

Next I will try the general problem (Et+Ez) again.

Eigensolver 
from dolfin import *
from mshr import *
import numpy as np


"""
# Define the forms - generates an generalized eigenproblem of the form
# [S]{h} = k_o^2[T]{h}
"""

def curl_t(w):
    return Dx(w[1], 0) - Dx(w[0], 1)



parameters["linear_algebra_backend"] = "PETSc"

def solver_conv(n_elements, a_r, b_r, electric=True):
    vector_order = 2
    nx = int(np.sqrt(n_elements))

    mesh = RectangleMesh(Point(0, 0), Point(a_r, b_r), nx, nx)
    dx = Measure('dx', domain=mesh)

    vector_element = FiniteElement("N1curl", mesh.ufl_cell(), vector_order)

    V = FunctionSpace(mesh, vector_element)

    # Define the test and trial functions
    u = TrialFunction(V)
    v = TestFunction(V)

    outer_boundary = DirichletBC(V, Constant((0.0, 0.0)), DomainBoundary())
    bcs = [outer_boundary]

    s = inner(curl_t(v), curl_t(u))*dx
    t = inner(v, u)*dx
    neigs = 40


    # Assemble the stiffness matrix (S) and mass matrix (T)
    S = PETScMatrix()
    T = PETScMatrix()
    assemble(s, tensor=S)
    assemble(t, tensor=T)

    if electric:
        outer_boundary.apply(S)
        outer_boundary.apply(T)
        [bc.zero(T) for bc in bcs]

    m = 1
    n= 0

    kc_squared = (m*pi/a_r)**2+(n*pi/b_r)**2

    # Solve the eigensystem
    solver = SLEPcEigenSolver(S, T)
    solver.parameters["spectrum"] = "target magnitude"
    solver.parameters["spectral_shift"] = kc_squared
    solver.parameters["spectral_transform"] = "shift-and-invert"
    solver.parameters["solver"] = "krylov-schur"

    solver.solve(neigs)

    return solver,V
18

Looks pretty good to me.

You should converge at some rate s such that \Vert u - u_h \Vert_{L_2} = \varepsilon \leq C h^s. You can find this rate of convergence by computing: \frac{\ln \varepsilon_c / \varepsilon_f}{\ln h_c / h_f} = s where the subscripts f and c denote a “fine” and a “coarse” representation, respectively.

According to the theory, for a linear problem you should converge at a rate s = p + 1 in the L_2 norm. Here’s a quick example with the Poisson equation which is hopefully easily adaptible:

from dolfin import *
from dolfin_dg import *
import numpy as np

p_order = 1
run_count = 0
ele_ns = [4, 8, 16, 32]
errorl2 = np.zeros(len(ele_ns))
errorh1 = np.zeros(len(ele_ns))
hsizes = np.zeros(len(ele_ns))

for j, ele_n in enumerate(ele_ns):
    mesh = UnitSquareMesh(ele_n, ele_n)

    V = FunctionSpace(mesh, 'CG', p_order)
    u, v = TrialFunction(V), TestFunction(V)

    gD = Expression('sin(pi*x[0])*sin(pi*x[1])',
                    domain=mesh, degree=p_order+1)

    f = -div(grad(gD))

    a = dot(grad(u), grad(v)) * dx
    L = f * v * dx

    bcs = [DirichletBC(V, gD, "on_boundary")]
    u = Function(V)
    solve(a == L, u, bcs)

    errorl2[j] = errornorm(gD, u, norm_type='l2', degree_rise=3)
    errorh1[j] = errornorm(gD, u, norm_type='h1', degree_rise=3)
    hsizes[j] = mesh.hmax()


if MPI.rank(mesh.mpi_comm()) == 0:
    l2rates = np.log(errorl2[:-1] / errorl2[1:]) / np.log(hsizes[:-1] / hsizes[1:])
    h1rates = np.log(errorh1[:-1] / errorh1[1:]) / np.log(hsizes[:-1] / hsizes[1:])
    print("L2 convergence rates: " + str(l2rates))
    print("H1 convergence rates: " + str(h1rates))

For the Hermitian eigenvalue problem I can’t recall exactly from memory what the rates should be. I think the eigenvalue itself converges at a rate of s_\lambda = 2p in the L_1 norm and the modes at a rate of s_\text{mode}=p in the L_2 norm.

Also if I recall correctly: in the non-Hermitian case, the rates of convergence are slower.

3 Likes
19

@nate Great, that helps a lot for understanding! I went back to my original problem, comparing with analytical solutions I get the expected modes/eigenvalues. However there is still plenty of non-physical eigenvalues around, what makes a convergence analysis tricky.

I went back to your original comment about the involution term, reading in literature that could definitely be the problem ( see for example: https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19950021305.pdf ).
When I plot div(E_t) I get scalar fields unequal to zero for the non-physical modes. Checking the derivation of the variational problem a zero divergence was never invoked.
I see that for curl-curl equations a zero divergence field can be achieved using Hodge-Laplacians (see https://fenicsproject.org/qa/10183/impose-div-a-0-for-curl-curl/, or Lagrangian multipliers to penalize for non-zero divergence.

Or should I just sacrifice computational efficiency and implement a div check on the results and disregard solutions with a non-zero divergence?

20

Thanks for this nice example, Nate. I like f = -div(grad(gD)) as I didn’t know something this straightforward was even possible. Would it be as easy to do a time derivative in a time-dependent problem (say the heat equation)?