Maxwell eigenvalue problem in 3D using Nédélec elements

dolfinx
1

Hi! I am trying to compute Maxwell eigenvalues in a cubical cavity [0, \pi]^3 using Nédélecs elements of the first kind, similar to this old demo which was recreated here for fenicsx/dolfinx.

At the end of this post I have attached my script where I have tried to adapt the code to do 3D instead of 2D. It seems to run fine, but I get incorrect eigenvalues (or at least I think I do!) I believe the 12 exact smallest nonzero eigenvalues should be:
Exact = [ 1.0 1.0 1.0 2.0 2.0 2.0 3.0 4.0 4.0 4.0 5.0 5.0]
but from my code I get:
Nédélec = [ 0.0 0.0 -0.0 2.0 2.0 2.0 3.0 3.0 5.0 5.0 5.0 5.0].
Am I expecting incorrect exact eigenvalues, or is my code incorrect somewhere? As far as I know, the extra dimension simply adds (due to the extra Fourier mode) another integer squared under the square root sign.

Code:

from slepc4py import SLEPc
from ufl import dx, curl, inner, TrialFunction, TestFunction
import numpy as np
from dolfinx.fem import (dirichletbc, Function, functionspace, form,
                         locate_dofs_topological)
from mpi4py import MPI
from dolfinx.fem.petsc import assemble_matrix
from petsc4py import PETSc
from dolfinx.mesh import (create_box, locate_entities_boundary,
                          CellType, GhostMode, exterior_facet_indices)
import sys


def par_print(string):
    if comm.rank == 0:
        print(string)
        sys.stdout.flush()

def eigenvalues(n_eigs, shift, V, bcs):
    # Define problem
    u, v = TrialFunction(V), TestFunction(V)
    a = form(inner(curl(u), curl(v)) * dx)
    b = form(inner(u, v) * dx)

    # Assemble matrices
    A = assemble_matrix(a, bcs, diagonal=1.0)
    A.assemble()
    # Zero rows of boundary DOFs of B. See [1]
    B = assemble_matrix(b, bcs, diagonal=0.0)
    B.assemble()

    # Create SLEPc Eigenvalue solver
    eps = SLEPc.EPS().create(comm)
    eps.setOperators(A, B)
    eps.setType(SLEPc.EPS.Type.KRYLOVSCHUR)
    eps.setProblemType(SLEPc.EPS.ProblemType.GHEP)
    eps.setWhichEigenpairs(eps.Which.TARGET_MAGNITUDE)
    eps.setTarget(shift)

    st = eps.getST()
    st.setType(SLEPc.ST.Type.SINVERT)
    st.setShift(shift)

    eps.setDimensions(n_eigs, PETSc.DECIDE, PETSc.DECIDE)
    eps.setFromOptions()
    eps.solve()

    its = eps.getIterationNumber()
    eps_type = eps.getType()
    n_ev, n_cv, mpd = eps.getDimensions()
    tol, max_it = eps.getTolerances()
    n_conv = eps.getConverged()

    par_print(f"Number of iterations: {its}")
    par_print(f"Solution method: {eps_type}")
    par_print(f"Number of requested eigenvalues: {n_ev}")
    par_print(f"Stopping condition: tol={tol}, maxit={max_it}")
    par_print(f"Number of converged eigenpairs: {n_conv}")

    computed_eigenvalues = []
    for i in range(min(n_conv, n_eigs)):
        lmbda = eps.getEigenvalue(i)
        computed_eigenvalues.append(np.round(np.real(lmbda), 1))

    eps.destroy()
    return np.sort(computed_eigenvalues)


def boundary(x):
    "Boundary marker"
    return boundary_lr(x) | boundary_tb(x) | boundary_fb(x)


def boundary_lr(x):
    "Left and right boundary marker"
    return np.isclose(x[0], 0.0) | np.isclose(x[0], np.pi)

def boundary_fb(x):
    "Front and back boundary marker"
    return np.isclose(x[1], 0.0) | np.isclose(x[1], np.pi)

def boundary_tb(x):
    "Top and bottom boundary marker"
    return np.isclose(x[2], 0.0) | np.isclose(x[2], np.pi)


def print_eigenvalues(mesh):
    # Nédélec
    V_n = functionspace(mesh, ("N1curl", 1))
    # Set boundary DOFs to 0 (u x n = 0 on \partial \Omega).

    # Version 1 utilizing "boundary" function
    # ud_n = Function(V_n)
    # f_dim = mesh.topology.dim - 1
    # boundary_facets = locate_entities_boundary(mesh, f_dim, boundary)
    # boundary_dofs_n = locate_dofs_topological(
    #     V_n, f_dim, boundary_facets)
    # bcs_nedelec = [dirichletbc(ud_n, boundary_dofs_n)]

    # Version 2 utilizing exterior_facet_indices
    bc_facets = exterior_facet_indices(mesh.topology)
    bc_dofs = locate_dofs_topological(V_n, mesh.topology.dim - 1, bc_facets)
    u_bc = Function(V_n)
    with u_bc.x.petsc_vec.localForm() as loc:
        loc.set(0)
    bcs_nedelec = [dirichletbc(u_bc, bc_dofs)]


    # Solve Maxwell eigenvalue problem
    eigenvalues_nedelec = eigenvalues(n_eigs, shift, V_n, bcs_nedelec)

    # Print results
    np.set_printoptions(formatter={'float': '{:5.1f}'.format})
    eigenvalues_exact = np.sort(np.array([float(m**2 + n**2 + k**2)
                                          for m in range(6)
                                          for n in range(6)
                                          for k in range(6)]))[1:13]
    par_print(f"Exact    = {eigenvalues_exact}")
    par_print(f"Nédélec  = {eigenvalues_nedelec}")

# Number of elements in each direction
n = 20
# Number of eigenvalues to compute
n_eigs = 12
# Find eigenvalues near
shift = 2.5
# Domain corners
corners = ((0.0, 0.0, 0.0), (np.pi, np.pi, np.pi))

comm = MPI.COMM_WORLD

par_print("Right diagonal mesh:")
mesh = create_box(
    comm,
    corners, (n, n, n),
    CellType.tetrahedron,
    ghost_mode=GhostMode.none)
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
print_eigenvalues(mesh)
2

I had a look through some texts other than the one I was using and actually the true eigenvalues are found with the fenicsx code. It is my interpretation of the reference which is wrong. A careful reading of Cheng “Field and wave electromagnetics” reveals that we should expect [2.0 2.0 2.0 3.0 3.0 5.0 5.0 5.0 5.0 5.0], since in our equation for the eigenvalues a^2 + b^2 + c^2, no two integers are allowed to be zero at the same time. Also multiplicity is more complicated as it has to do with the vectors spanning the plane perpendicular to the normal vector (a,b,c), this is why 3 appears twice though there is only one way to get 3=1+1+1. See https://arxiv.org/pdf/1306.5354. I leave this for posterity in case someone else stumbles across this problem

2 Likes