Reproducing Stable and unstable finite elements for the Maxwell eigenvalue problem in Fenicsx

Linear Algebra
1

I am trying to reproduce this tutorial about stable and unstable finite elements for Maxwell’s equations in Fenicsx. However, I am struggling with the set up of the SLEPc solver. So far, I have the following MWE, where only the Lagrangian elements are working yet:

import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from slepc4py import SLEPc
from dolfinx import mesh
from dolfinx.fem import FunctionSpace, VectorFunctionSpace, form, dirichletbc, locate_dofs_geometrical, locate_dofs_topological
from dolfinx.fem.petsc import assemble_matrix
from ufl import FiniteElement, VectorElement, TestFunction, TrialFunction, curl, dot, dx, grad, inner


def eigenvalues(V, bcs):

    u = TrialFunction(V)
    v = TestFunction(V)
    a = form(inner(curl(u), curl(v))*dx)
    b = form(inner(u, v)*dx)

    A = assemble_matrix(a, bcs)
    A.assemble()
    B = assemble_matrix(b, bcs)
    B.assemble()
    # [bc.zero(B) for bc in bcs] # in the fenics tutorial they zero out the corresponding lines, why?


    solver = SLEPc.EPS()
    solver.create()
    solver.setOperators(A, B)
    solver.setProblemType(SLEPc.EPS.ProblemType.GHEP)
    solver.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_MAGNITUDE)
    solver.setTarget(5.5)
    neigs = 12
    solver.setDimensions(nev=neigs, ncv=2*neigs)
    solver.setFromOptions()

    solver.solve()
    print(solver.getConverged())
    # Return the computed eigenvalues in a sorted array
    computed_eigenvalues = []
    for i in range(min(neigs, solver.getConverged())):
        r = solver.getEigenvalue(i)  # ignore the imaginary part
        computed_eigenvalues.append(r)
    return np.sort(np.array(computed_eigenvalues))


ncells = 40
corners = [[0.0, 0.0], [np.pi, np.pi]]
msh = mesh.create_rectangle(MPI.COMM_WORLD, corners, [ncells, ncells])

V = VectorFunctionSpace(msh, ("CG", 1))
# V = FunctionSpace(msh, ("N1curl", 1))

def top(x):
    return np.logical_and(
        np.isclose(x[1], np.pi), x[0] < np.pi
    )

def bottom(x):
    return np.logical_and(
        np.isclose(x[1], 0.0), x[0] > 0
    )

def right(x):
    return np.logical_and(
        np.isclose(x[0], np.pi), x[1] < np.pi
    )

def left(x):
    return np.logical_and(
        np.isclose(x[0], 0.0), x[1] > 0
    )

bd_markers_y = [top, bottom]
bd_markers_x = [left, right]
bcs = []

for bd_marker in bd_markers_y:
    bd_facets = mesh.locate_entities_boundary(msh, msh.topology.dim-1, bd_marker)
    bd_dofs_x = locate_dofs_topological(V.sub(0), msh.topology.dim-1, bd_facets)
    bc = dirichletbc(ScalarType(0), bd_dofs_x, V.sub(0))
    bcs.append(bc)

for bd_marker in bd_markers_x:
    bd_facets = mesh.locate_entities_boundary(msh, msh.topology.dim - 1, bd_marker)
    bd_dofs_y = locate_dofs_topological(V.sub(1), msh.topology.dim - 1, bd_facets)
    bc = dirichletbc(ScalarType(0), bd_dofs_y, V.sub(1))
    bcs.append(bc)

ev = eigenvalues(V, bcs)
print(ev)

More precisely, I have two questions

  1. In the tutorial, the lines corresponding to the BCs in the operator B are zeroed out. How can this be done in Fenicsx?
  2. When I comment out the lines
solver.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_MAGNITUDE)
solver.setTarget(5.5)

I get the 12 eigenvalues with the largest magnitude as expected. However, when I include the lines, no eigenvalue converges. What could be the cause of this?

2

I can only answer this point of your question.
You can set the keyword argument diagonal=0 in assemble_matrix, dolfinx.fem.petsc — DOLFINx 0.5.1 documentation
i.e.

B = assemble_matrix(b, bcs, diagonal=0)

As you have not supplied a reproducible example for the other parts of the code, I cannot really give you any further guidance (i.e. I cannot run your code, as Im missing imports and definitions of variables).

3

As far as I know, setWhichEigenpairs and setTarget are useful when using spectral transformation in SLEPc, which you are not using here. To do that, add the following lines to your code:

solver = SLEPc.EPS()
solver.create()
solver.setOperators(A, B)

# Add this #
solver.setType(SLEPc.EPS.Type.KRYLOVSCHUR)
st = solver.getST()
st.setType(SLEPc.ST.Type.SINVERT)
############

solver.setProblemType(SLEPc.EPS.ProblemType.GHEP)
solver.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_MAGNITUDE)
solver.setTarget(5.5)
neigs = 12
solver.setDimensions(nev=neigs, ncv=2*neigs)
solver.setFromOptions()

Also, this post may be useful:

1 Like
4

I updated this demo here: GitHub - jpdean/maxwell_eigenvalue

3 Likes
5

@jpdean Quick follow-up question: In the old version, there was the following line

[bc.zero(B) for bc in bcs]

The explanation for it was

We zero out the rows of B
corresponding to constrained boundary degrees of freedom, so as not to introduce spurious eigenpairs with nonzero boundary DOFs.

In your implementation, this translates to

# Zero rows of boundary DOFs of B. See [1]
B = assemble_matrix(b, bcs, diagonal=0.0)

Could you quickly explain, what this does exactly and why it is equivalent?

6

When assemble_matrix is called on a form where the trial and test spaces are the same, the rows and columns of the degrees of freedom constrained by Dirichlet boundary conditions are zeroed out and, by default, a 1 is placed on the diagonal. This ensures that the constrained degrees of freedom have the correct value when the right hand side vector is modified appropriately (i.e. by apply_lifting and set_bc). Here however, we don’t want ones on the diagonal, so we pass the optional argument diagonal = 0.0 when calling assemble_matrix.

4 Likes