Maxwell solver example does not work after 0.5 update

dolfinx
1

An old Maxwell’s equation solver example (Stable and unstable finite elements for the Maxwell eigenvalue problem — DOLFIN documentation) ported to DolfinX (Waveguide number of eigenmodes - #8 by MatsievskiySV) used to work fine prior to v0.5, but now it does not produce physical solutions. The only difference is substitution of compute_boundary_facets with locate_entities_boundary. What is wrong here?

Full code:

from dolfinx import (io, fem, mesh)
import ufl
import numpy as np
from mpi4py import MPI
from slepc4py import SLEPc
from petsc4py.PETSc import ScalarType


msh = mesh.create_rectangle(comm=MPI.COMM_WORLD,
                            points=((0, 0), (2, 1)),
                            n=(40, 40),
                            cell_type=mesh.CellType.triangle)

msh.topology.create_connectivity(msh.topology.dim-1, msh.topology.dim)

N1curl = ufl.FiniteElement("Nedelec 1st kind H(curl)", msh.ufl_cell(), 1)
V = fem.FunctionSpace(msh, N1curl)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

a = ufl.inner(ufl.curl(u), ufl.curl(v)) * ufl.dx
b = ufl.inner(u, v) * ufl.dx

bc_facets = mesh.locate_entities_boundary(msh,
                                          msh.topology.dim - 1,
                                          lambda x: np.full(x.shape[1], True, dtype=bool))
bc_dofs = fem.locate_dofs_topological(V, msh.topology.dim-1, bc_facets)

u_bc = fem.Function(V)
u_bc.x.array[:] = ScalarType(0)
bc = fem.dirichletbc(u_bc, bc_dofs)

A = fem.petsc.create_matrix(fem.form(a))
fem.petsc.assemble_matrix(A, fem.form(a), bcs=[bc])
A.assemble()
B = fem.petsc.create_matrix(fem.form(b))
fem.petsc.assemble_matrix(B, fem.form(b), bcs=[bc])
B.assemble()

eps = SLEPc.EPS().create(MPI.COMM_WORLD)
eps.setOperators(A, B)
eps.setType(SLEPc.EPS.Type.KRYLOVSCHUR)
st = eps.getST()
st.setType(SLEPc.ST.Type.SINVERT)
eps.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_MAGNITUDE)
eps.setTarget(12)
eps.setDimensions(30)
eps.solve()

vals = [(i, eps.getEigenvalue(i)) for i in range(eps.getConverged())
        if not np.isclose(eps.getEigenvalue(i), 0)
        and not np.isclose(eps.getEigenvalue(i), 1)]
vals.sort(key=lambda x: x[1].real)

j = 0
E = fem.Function(V)
with io.XDMFFile(MPI.COMM_WORLD, "out.xdmf", "w") as xdmf:
    xdmf.write_mesh(msh)
    for i, _ in vals:
        print(f"eigenvalue: {eps.getEigenpair(i, E.vector).real:.12f}")
        E.name = f"E-{j+1:03d}-{eps.getEigenpair(i, E.vector).real:.4f}"
        xdmf.write_function(E)
        j += 1

OS: Debian bookworm
Built from main branches as of 24-08-2022:
BASIX commit: Update UFL wrapper (#575) · FEniCS/basix@1fdeeb4 · GitHub
DOLFIN commit: correct gdim/tdim (#2307) · FEniCS/dolfinx@e3d7b51 · GitHub
UFL commit: Sanitise unicode input (#111) · FEniCS/ufl@50d2a48 · GitHub
FFCX commit: Use Basix element functionality over UFL (#515) · FEniCS/ffcx@4bc7476 · GitHub

2

Could you add what eigenvalues you expect to get printed?
When I run the code below in v0.4.1 i get:

eigenvalue: 2.467453845009
eigenvalue: 9.860312654523
eigenvalue: 9.870441650323
eigenvalue: 12.331992693807
eigenvalue: 19.747631128070
eigenvalue: 22.210852229644
eigenvalue: 32.109120054580
eigenvalue: 39.330212890661
eigenvalue: 39.491630029626
eigenvalue: 41.814442297205
eigenvalue: 49.267973337875
eigenvalue: 49.419419448067
eigenvalue: 61.692726728336
eigenvalue: 61.716767616632
eigenvalue: 71.681859411201
eigenvalue: 79.092163297278
eigenvalue: 88.079821259095

while if i run it with v0.5.0 I get

eigenvalue: 2.467453845010
eigenvalue: 9.860312654523
eigenvalue: 9.870441650323
eigenvalue: 12.331992693807
eigenvalue: 19.747631128070
eigenvalue: 22.210852229644
eigenvalue: 32.109120054579
eigenvalue: 39.330212890658
eigenvalue: 39.491630029626
eigenvalue: 41.814442297208
eigenvalue: 49.267973337872
eigenvalue: 49.419419448067
eigenvalue: 61.692726728332
eigenvalue: 61.716767616630
eigenvalue: 71.681859411199
eigenvalue: 79.092163297276

where the only difference (up to the 6-7 digit is that v0.4.1 finds one more eigenvalue).

Note that you could use

bc_facets = mesh.exterior_facet_indices(msh.topology)

in v0.5.0, producing the same result as v0.5.0 above.

Please note that you should never use XDMFFile for post-processing Nedelec-kind elements (as it interpolates into CG-1, which is not well defined).

You should use VTXWriter after interpolating the solution into DG-1.

3

This is weird. After rebuilding from the same commits everything is back to normal.

This is very handy, thanks! I wish your changelog was up-to-date, though.

4

Its quite hard for me to keep an extensive changelog between releases.
A summary can be found at