Periodic BCs in Dolfinx with eigensolver

Hello all,
I am presently attempting to solve the following weak form in Dolfinx (which I have successfully solved using legacy Dolfin, which is discussed here):
F=\left (\Gamma\hat{\mathbf{z}}\times\mathbf{E} + \nabla\times\mathbf{E}\right )^*\cdot\left (\Gamma\hat{\mathbf{z}}\times\mathbf{E}+\nabla\times\mathbf{E}\right ) - k_0^2\mathbf{E}^*\cdot\mathbf{E},
where the periodic boundary is enforced between z=0 and z=l. The term \Gamma is the complex valued Floquet propagation factor (for lossless structures, we assume to take imaginary values only).

I have a couple of problems -

  1. In line 132 of the code (where I define the MPC constraint), there is a problem with the variable types, but after looking at the man page for “create_periodic_constraint_topological” and looking at examples, something is still not right. Perhaps my boundary relation function output is not correctly presented.
  2. What is the best way to set up the resulting eigenvalue problem?
# Generate unit cell of helical structure
import sys
import numpy as np
from mpi4py import MPI as nMPI
from petsc4py import PETSc
import slepc4py as slepc
import meshio
import ufl
from dolfinx.mesh import locate_entities_boundary, meshtags
from dolfinx.io import gmshio
from dolfinx import fem
from dolfinx import io
import dolfinx_mpc
import gmsh

l = 6.0 # Cell length
r = 6.0 # Cell radius
rh = 2.0 # Helix radius
rw = 0.2 # wire radius
npts = 100 # number of points of spline curve
eps = 1.0e-5

Gamma = 0 + 6.283j

comm = nMPI.COMM_WORLD
mpiRank = comm.rank
modelRank = 0 # process where GMSG runs

if mpiRank == modelRank:

    gmsh.initialize()
    gmsh.option.setNumber('General.Terminal', 1)
    gmsh.model.add("Helix Cell")

    c = gmsh.model.occ.addCylinder(0, 0, 0, 0, 0, l, r, 1, 2*np.pi)

    dt = 2 * np.pi / npts
    p = []
    for n in range(npts+3):
        theta = 2 * np.pi * (n-1) / npts
        kk = gmsh.model.occ.addCircle(rh * np.cos(theta), rh * np.sin(theta), (n -1) * l / npts, rw)
        p.append(gmsh.model.occ.addCurveLoop([kk], tag=-1))
    print(p)
    s = gmsh.model.occ.addThruSections(p, tag=-1)
    print(s)
    v = gmsh.model.occ.cut([(3,c)], s, tag=100, removeObject=True, removeTool=True)
    print("v={}".format(v))
    for n in range(npts+3):
        gmsh.model.occ.remove([(1,p[n])], recursive=True)

    gmsh.model.occ.synchronize()

    bb = gmsh.model.getEntities(dim=3)
    pt = gmsh.model.getBoundary(bb, combined=True, oriented=False, recursive=False)
    print(pt)
# Boundary conditions
    Bot = []
    Top = []
    PEC = []
    for bb in pt:
        CoM = gmsh.model.occ.getCenterOfMass(bb[0], bb[1])
        print(CoM)
        if(np.abs(CoM[2]) < eps):
            Bot.append(bb[1])
        elif(np.abs(CoM[2]-l) < eps):
            Top.append(bb[1])
        else:
            PEC.append(bb[1])

    print(Bot, Top, PEC)
    gmsh.model.mesh.setPeriodic(2, Top, Bot, [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, l, 0, 0, 0, 1])

    gmsh.model.addPhysicalGroup(3, [100], 1) # Helix volume
    gmsh.model.setPhysicalName(3, 1, "HelixVolume")
    gmsh.model.addPhysicalGroup(2, Top, 1) # Top BC (slave periodic BC)
    gmsh.model.addPhysicalGroup(2, Bot, 2) # Bottom BC (master periodic BC)
    gmsh.model.addPhysicalGroup(2, PEC, 3) # PEC boundary
    gmsh.model.setPhysicalName(2, 1, "TopPBC")
    gmsh.model.setPhysicalName(2, 2, "BotPBC")
    gmsh.model.setPhysicalName(2, 3, "PEC")

    gmsh.model.mesh.setSize([(0,108),(0,109)], 0.05)
    gmsh.option.setNumber('Mesh.MeshSizeMin', 0.001)
    gmsh.option.setNumber('Mesh.MeshSizeMax', 1.2)
    gmsh.option.setNumber('Mesh.Algorithm', 6) #1=Mesh Adapt, 2=Auto, 3=Initial mesh only, 5=Delaunay, 6=Frontal-Delaunay
    gmsh.option.setNumber('Mesh.MshFileVersion', 2.2)
    gmsh.option.setNumber('Mesh.Format', 1)
    gmsh.option.setNumber('Mesh.MinimumCirclePoints', 50)
    gmsh.option.setNumber('Mesh.CharacteristicLengthFromCurvature', 1)
    gmsh.model.mesh.generate(3)
    gmsh.fltk.run()

# mesh = 3d topology
# ct = volume tags
# fm = surface tags
mesh, ct, fm = gmshio.model_to_mesh(gmsh.model, comm, modelRank, gdim=3)
if modelRank == mpiRank:
    gmsh.finalize()

mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
elem = ufl.FiniteElement('Nedelec 1st kind H(curl)', mesh.ufl_cell(), degree=2)
V = fem.FunctionSpace(mesh, elem)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# Print out BC tags
with io.XDMFFile(mesh.comm, "BCs.xdmf", "w") as xx:
    xx.write_mesh(mesh)
    xx.write_meshtags(fm)

n = ufl.FacetNormal(mesh)

# Periodic BC mapping
def PeriodicMapping(x):
    out = np.zeros(x.shape)
    out[0] = x[0]
    out[1] = x[1]
    out[2] = x[2] + l
    return out

# Dirichlet BCs
facets = fm.find(3) # PEC facets
ubc = fem.Function(V)
ubc.x.set(0+0j)
dofs = fem.locate_dofs_topological(V, mesh.topology.dim-1, facets)
bc = fem.dirichletbc(ubc, dofs)

# Create multipoint boundary.  Master boundary has tag 2, 
# slave has tag 1 - set on periodic constraint.
# (see GMSH physical group definitions above)
mpc = dolfinx_mpc.MultiPointConstraint(V)
mpc.create_periodic_constraint_topological(V, fm, 1, PeriodicMapping, bc, 1.0)
mpc.finalize()

# Set up forms for eigenvalue problem Ax -k0^2*Bx = 0
az = ufl.asvector((0, 0, 1)) # Z-directed unit vector
a = ufl.inner(Gamma * ufl.cross(az, u) + ufl.curl(u), Gamma * ufl.cross(az, v) + ufl.curl(v)) * dx
b = ufl.inner(u, v) * dx

A = dolfinx_mpc.assemble_matrix(a, mpc, bc, 1.0)
A.assemble()
B = dolfinx_mpc.assemble_matrix(b, mpc, bc, 1.0)
B.assemble()
eigs = slepc.EPS().create(mesh.comm)
eigs.setOperators(A, B)
eigs.setProblemType(slepc.EPS.ProblemType.GNHEP) # General non-hermitian
eigs.setTolerances(1.0e-8)
eigs.setType(slepc.EPS.Type.KRYLOVSCHUR)
st = eigs.getST() # Invoke spetcral transformation
st.setType(slepc.ST.Type.SINVERT)
eigs.setWhichEigenpairs(slepc.EPS.Which.TARGET_REAL) # target real part of eigenvalue
eigs.setTarget(-1.0)
eigs.setDimensions(nev=5) # Number of eigenvalues
eigs.solve()
eigs.view()
eigs.errorView()

sys.exit(0)

You are using Nedelec-elements, which has not been tested extensively with DOLFINX_MPC, as the functionals are not simple point evaluations (but integral moments), which makes setting up the relations quite tricky (as one has to consider, what does a periodic condition mean in general for Nedelec elements.

The periodicity would be well-defined if the facets on either end are lining up (i.e. a facet on the left side maps perfectly to a facet on the right side).

That’s food for thought. Would it matter that the GMSH mesh generator has generated a periodic mesh? (This was done.) The mesh on both faces should be equal. I used Nedelec elements in the legacy Dolfin model with no problem with a GMSH periodic mesh.

How can I explicitly treat the facets on left and right sides (as these are available from the GMSH physical surface data)?

In old Dolfin Nedelec elements were based on point evaluations, causing sub-optimal interpolation results, ref: fenics-project / FIAT / issues / #25 - Nedelec DOFs give suboptimal interpolation operator — Bitbucket

The implementation of periodic conditions in legacy dolfin is also very different from the one done in MPC.
In legacy DOLFIN the dofmap is modified, while in DOLFINx_mpc, the matrices are modified at assembly.

I would have to think about if the periodic conditions for perfectly matched meshes works as of now.
I would suggest making an issue Sign in to GitHub · GitHub so I dont forget it.

Thanks. Opened issue in Github
Cheers!

I think I may try to formulate the problem using Lagrange interpolation to start off. I changed to a vector lagrange interpolation and MPC seemed to work. Question: I will need to invoke a Dirichlet condition on the perfect electric conducting surfaces. What is the best way to do \mathbf{n}\times\mathbf{E}=0 along a surface (where n is the normal unit vector along the Dirichlet boundary)? The boundaries are not generally aligned with the vector component directions.

UPDATE: I will try to use a weak form of the Dirichlet BCs, which I use in this version of the code.

Hi,

I also need to use periodic boundaries for Nedelec elements. @dokken Is there any update in the development branch?

@BillS Could you share what came out of your experiment with vector Lagrange elements?

Thanks for your help.

Hello,
I have not been able to get this to work on my installation yet. Unfortunately, I have had to dedicate time to other things. I am not sure when I will be able to look at this problem again, but hopefully I can before the end of the year.