Using Sherman Morrison formula for generating rank-1 boundary operator with Fenicsx and PETSc

Linear Algebra
1

Good day all,

I finally have got around to looking at modeling a generalized waveguide port boundary condition. Starting with the simplest case of a single mode, I am trying to implement the scattered port wave operator, which is a rank-1 projection operator, whose contribution to the full weak form is

B(\mathbf{E}) = \frac{jk_0\eta_0}{Z_0}\int\!\!\!\int\limits_{\partial\Omega}\mathbf{E}^*\times\mathbf{n}\cdot\mathbf{h}_0 dS \int\!\!\!\int\limits_{\partial\Omega}\mathbf{E}\times\mathbf{n}\cdot\mathbf{h}_0 \, dS

where Z_0 is the mode impedance of the TE10 mode and \mathbf{h}_0 is the transverse magnetic field for the TE10 mode.

Unless I am mistaken, the operator above can be represented by the outer product of the vector produced by

u=\int\!\!\!\int\limits_{\partial\Omega}\mathbf{v}\times\mathbf{n}\cdot\mathbf{h}_0 \, dS

Of course, constructing the dense outer product makes no sense, so I am trying to implement the Sherman-Morrison method of solving the matrix system

(A+uu^T) x = b

where u is the vector generated by the integral form above on each of the waveguide ports.

Here is the code. The part of interest is at the end (where I do the vector manipulation).

# Example to test multi-mode waveguide absorbing boundary condition
# for a simple rectangular waveguide
import sys
import numpy as np
from mpi4py import MPI as M
from petsc4py import PETSc as P
import gmsh
import ufl as U
import basix.ufl
from dolfinx.mesh import locate_entities_boundary, meshtags
from dolfinx.io import gmshio
from dolfinx.fem import petsc
from dolfinx import fem
from dolfinx import io
import time

eta0 = 377.0 # free space wave impedance

# The TE modes in rectangular WG for excitation and mode vectors in absorbing BC
class ModesTE:
    def __init__(self, beta, k, a, b, m, n):
        self.beta = beta
        self.k = k
        self.a = a
        self.b = b
        self.m = m
        self.n = n
    def eval_H(self, x):
        if((self.m != 0) and (self.n != 0)):
            B = np.sqrt(4.0 / (self.a * self.b))
        elif self.m == 0:
            B = np.sqrt(2.0 / (self.a * self.b))
        else:
            B = np.sqrt(2.0 / (self.a * self.b))
        hx = B * np.sin((np.pi * self.m / self.a) * x[0]) * np.cos((np.pi * self.n / self.b) * x[1])
        hy = B * np.cos((np.pi * self.m / self.a) * x[0]) * np.sin((np.pi * self.n / self.b) * x[1])
        hz = np.full_like(hx, 0.0+0j, dtype=np.complex128)
        return(hx, hy, hz)

eps = 1.0e-4
a = 22.86 #WG width
b = 10.16 #WG height
l = 50.8 # WG length
f = 10.5 #frequency of operation GHz
c = 300.0 #speed of light
lm = 1.60

comm = M.COMM_WORLD
mpiRank = comm.rank
modelRank = 0 # where GMSH runs
k0 = 2.0 * np.pi * f / c #wavenumber wrt frequency

if mpiRank == modelRank:
    gmsh.initialize()
    gmsh.option.setNumber('General.Terminal', 1)
    gmsh.model.add('Waveguide')
    gmsh.model.occ.addBox(0, 0, 0, a, b, l, 1)
    gmsh.model.occ.synchronize()
    pt = gmsh.model.getEntities(dim=0)
    gmsh.model.mesh.setSize(pt, lm)

    gmsh.option.setNumber('Mesh.MeshSizeMin', 0.01)
    gmsh.option.setNumber('Mesh.MeshSizeMax', 2.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.MinimumCirclePoints', 30)
    gmsh.option.setNumber('Mesh.CharacteristicLengthFromCurvature', 1)
    gmsh.option.setNumber('Mesh.Tetrahedra', 0)
    OUT = []
    IN = []
    PEC = []
    bb = gmsh.model.getEntities(dim=3)
    pt = gmsh.model.getBoundary(bb, combined=True, oriented=False, recursive=False)
    for bnd in pt:
        CoM = gmsh.model.occ.getCenterOfMass(bnd[0], bnd[1])
        if np.abs(CoM[2]) < eps:
            IN.append(bnd[1])
        elif np.abs(CoM[2] - l) < eps:
            OUT.append(bnd[1])
        else:
            PEC.append(bnd[1])

    gmsh.model.addPhysicalGroup(3, [1], 1) #WG
    gmsh.model.setPhysicalName(3, 1, "Waveguide")
    gmsh.model.addPhysicalGroup(2, PEC, 1)
    gmsh.model.setPhysicalName(2, 1, "PEC")
    gmsh.model.addPhysicalGroup(2, IN, 2)
    gmsh.model.setPhysicalName(2, 2, "INPORT")
    gmsh.model.addPhysicalGroup(2, OUT, 3)
    gmsh.model.setPhysicalName(2, 3, "OUTPORT")

    gmsh.model.mesh.generate(3)
#    gmsh.fltk.run()


mesh, ct, fm = gmshio.model_to_mesh(gmsh.model, comm, modelRank, gdim=3)
if mpiRank == modelRank:
    gmsh.finalize()

mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
elem = basix.ufl.element('Nedelec 1st kind H(curl)', mesh.basix_cell(), degree=2)
V = fem.functionspace(mesh, elem)
u = U.TrialFunction(V)
v = U.TestFunction(V)
normal = U.FacetNormal(mesh)
ds = U.Measure("ds", domain=mesh, subdomain_data = fm)
# Dirichlet BCs on waveguide walls only in-ports and out ports have Neumann condition.
facets = fm.find(1)
ubc = fem.Function(V)
ubc.interpolate(lambda x: np.array([0*x[0], 0*x[1], 0*x[2]], dtype=np.complex128))
dofs = fem.locate_dofs_topological(V, mesh.topology.dim-1, facets)
bc = fem.dirichletbc(ubc, dofs)
# Build rhs - Excitation at in port.
uh = fem.Function(V, dtype=np.complex128)
kz = np.sqrt(k0 * k0 - (np.pi / a)**2)
f = ModesTE(kz, k0, a, b, 1, 0)
uh.interpolate(f.eval_H)
finc = 2.0j * k0 * eta0 * U.inner(U.cross(normal, uh), v) * ds(2)
uh.name = "H_inc"
# Generate mode inner product for projection operator overi in and out port surfaces.
ua = fem.Function(V, dtype=np.complex128)
kz = np.sqrt(k0 * k0 - (np.pi / a)**2)
Zm = eta0 * k0 / kz
f = ModesTE(kz, k0, a, b, 1, 0) #Mode 1 only
ua.interpolate(f.eval_H)
# Generate linear forms for mode functions on ports
um1a = fem.form(U.inner(ua, U.cross(normal, v)) * ds(2)) #input plane
um1b = fem.form(U.inner(ua, U.cross(normal, v)) * ds(3)) #output plane

t_start = time.perf_counter()

# Weak form
L = (U.inner(U.curl(u), U.curl(v)) - k0 * k0 * U.inner(u, v)) * U.dx
Df = fem.form(L)
Rf = fem.form(finc)

num_dofs_local = V.dofmap.index_map.size_local * V.dofmap.index_map_bs
num_dofs_global = V.dofmap.index_map.size_global * V.dofmap.index_map_bs
print(f"Number of dofs (owned) by rank {mpiRank}: {num_dofs_local}")
if mpiRank == modelRank:
    print(f"Number of dofs global: {num_dofs_global}")

#Define PETSc vectors for mode projection operator of form UU^T
def VecCreate(linear_form, bilin_form, BCs):
    vec = petsc.create_vector(linear_form)
    with vec.localForm() as lb:
        lb.set(0.0)
    petsc.assemble_vector(vec, linear_form)
    petsc.apply_lifting(vec, [bilin_form], [[BCs]])
    vec.ghostUpdate(addv=P.InsertMode.ADD_VALUES, mode=P.ScatterMode.REVERSE)
    petsc.set_bc(vec, [BCs])
    vec.ghostUpdate(addv=P.InsertMode.INSERT, mode=P.ScatterMode.FORWARD)
    return vec

#A = petsc.create_matrix(Df)
#A.zeroEntries()
A = petsc.assemble_matrix(Df, bcs=[bc])
A.assemble()
b = VecCreate(Rf, Df, bc)
u1 = VecCreate(um1a, Df, bc)
u2 = VecCreate(um1b, Df, bc)

petsc_options = {
        "ksp_type": "preonly",
        "pc_type": "lu",
        "pc_factor_mat_solver": "mumps",
        "pc_factor_nonzeros_along_diagonal": 1e-8,
        "ksp_error_if_not_converged": False,
        "ksp_monitor": None }
opts = P.Options()
for key, value in petsc_options.items():
    opts[key] = value

solver = P.KSP().create(mesh.comm)
solver.setOperators(A)
solver.setFromOptions()
solver.setUp()
t_stop = time.perf_counter()

print("Solver setup time = {0}".format(t_stop - t_start))
sol = fem.Function(V)
sol.name="Geom Solution"
w1 = fem.Function(V)
w1.name="BC1"
w2 = fem.Function(V)
w2.name="BC2"

# Solve base problem A*sol=b 
solver.solve(b, sol.x.petsc_vec)
sol.x.scatter_forward()
t_solve = time.perf_counter()

print("Solve time = {0}".format(t_solve - t_stop))
# Generate A*w1=u1 for Sherman Morrison on in port.
solver.solve(u1, w1.x.petsc_vec)
w1.x.scatter_forward()
t_solve1 = time.perf_counter()
print("Solve1 time = {0}".format(t_solve1 - t_solve))
# Generate A*w2=u2 for Sherman morrison on out port
solver.solve(u2, w2.x.petsc_vec)
w2.x.scatter_forward()
t_solve2 = time.perf_counter()

print("Solve1 time = {0}".format(t_solve2 - t_solve1))

#Perform PETSc vector operations to construct Sherman Morrison solution
#                 D1             D2
# fsol= sol - --------- w1 - ----------w2
#               (1+C1)         (1+C2)
#
kz = np.sqrt(k0 * k0 - (np.pi / a)**2)
C1 = u1.tDot(w1.x.petsc_vec)
C2 = u2.tDot(w2.x.petsc_vec) 
D1 = u1.tDot(sol.x.petsc_vec)
D2 = u2.tDot(sol.x.petsc_vec)
print(C1, C2, D1, D2)
# Compute the composite solution with WG BC contributions
fsol = fem.Function(V)
fsol.name="Full Solution"
fsol.x.petsc_vec.set(0.0)
fsol.x.petsc_vec.axpy(1.0, sol.x.petsc_vec) # Solution without rad BCs
fsol.x.petsc_vec.axpy(-1.0j * kz * D1 / (1.0 + C1), w1.x.petsc_vec) # Input rad BC
fsol.x.petsc_vec.axpy(-1.0j * kz * D2 / (1.0 + C2), w2.x.petsc_vec) # Output rad BC

#Plot!
Vw = basix.ufl.element('DG', mesh.basix_cell(), 0, shape=(mesh.geometry.dim, ))
W = fem.functionspace(mesh, Vw)
Et = fem.Function(W)
Et.interpolate(fsol)
Et.name = "ElectricField"
with io.XDMFFile(mesh.comm, "Efield.xdmf", "w") as xx:
    xx.write_mesh(mesh)
    xx.write_function(Et)

sys.exit(0)

The code runs without problems, but there is an error that causes the code to return a solution that is not correct. I am not sure if I am treating the PETSc vector operations correctly (still being a bit of an newbie, in this regard). The real and imaginary parts of the solution should have equal magnitude and be 90 degrees out of phase. Using Paraview to observe the result shows that this is not the case. The result should be closely equivalent to the result of the code below (where I use the waveguide absorbing BC that I have been using up to now - that gives the corerct solution.)

# Example to test multi-mode waveguide absorbing boundary condition
# for a simple rectangular waveguide
import sys
import numpy as np
from mpi4py import MPI as M
from petsc4py import PETSc as P
import gmsh
import ufl as U
import basix.ufl
from dolfinx.mesh import locate_entities_boundary, meshtags
from dolfinx.io import gmshio
from dolfinx.fem import petsc
from dolfinx import fem
from dolfinx import io
import time

eta0 = 377.0 # free space wave impedance

class ModesTE:
    def __init__(self, beta, k, a, b, m, n):
        self.beta = beta
        self.k = k
        self.a = a
        self.b = b
        self.m = m
        self.n = n
    def eval_H(self, x):
        if((self.m != 0) and (self.n != 0)):
            B = np.sqrt(self.n * self.m * 4.0 / (self.a * self.b))
        elif self.m == 0:
            B = np.sqrt(self.n * 2.0 / (self.a * self.b))
        else:
            B = np.sqrt(self.m * 2.0 / (self.a * self.b))
        hx = B * np.sin((np.pi * self.m / self.a) * x[0]) * np.cos((np.pi * self.n / self.b) * x[1])
        hy = B * np.cos((np.pi * self.m / self.a) * x[0]) * np.sin((np.pi * self.n / self.b) * x[1])
        hz = np.full_like(hx, 0.0+0j, dtype=np.complex128)
        return(hx, hy, hz)

eps = 1.0e-4
a = 22.86 #WG width
b = 10.16 #WG height
l = 50.8 # WG length
f = 10.5 #frequency of operation GHz
c = 300.0 #speed of light
lm = 1.60

comm = M.COMM_WORLD
mpiRank = comm.rank
modelRank = 0 # where GMSH runs
k0 = 2.0 * np.pi * f / c #wavenumber wrt frequency

if mpiRank == modelRank:
    gmsh.initialize()
    gmsh.option.setNumber('General.Terminal', 1)
    gmsh.model.add('Waveguide')
    gmsh.model.occ.addBox(0, 0, 0, a, b, l, 1)
    gmsh.model.occ.synchronize()
    pt = gmsh.model.getEntities(dim=0)
    gmsh.model.mesh.setSize(pt, lm)

    gmsh.option.setNumber('Mesh.MeshSizeMin', 0.01)
    gmsh.option.setNumber('Mesh.MeshSizeMax', 2.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.MinimumCirclePoints', 30)
    gmsh.option.setNumber('Mesh.CharacteristicLengthFromCurvature', 1)
    gmsh.option.setNumber('Mesh.Tetrahedra', 0)
    OUT = []
    IN = []
    PEC = []
    bb = gmsh.model.getEntities(dim=3)
    pt = gmsh.model.getBoundary(bb, combined=True, oriented=False, recursive=False)
    for bnd in pt:
        CoM = gmsh.model.occ.getCenterOfMass(bnd[0], bnd[1])
        if np.abs(CoM[2]) < eps:
            IN.append(bnd[1])
        elif np.abs(CoM[2] - l) < eps:
            OUT.append(bnd[1])
        else:
            PEC.append(bnd[1])

    gmsh.model.addPhysicalGroup(3, [1], 1) #WG
    gmsh.model.setPhysicalName(3, 1, "Waveguide")
    gmsh.model.addPhysicalGroup(2, PEC, 1)
    gmsh.model.setPhysicalName(2, 1, "PEC")
    gmsh.model.addPhysicalGroup(2, IN, 2)
    gmsh.model.setPhysicalName(2, 2, "INPORT")
    gmsh.model.addPhysicalGroup(2, OUT, 3)
    gmsh.model.setPhysicalName(2, 3, "OUTPORT")

    gmsh.model.mesh.generate(3)
    gmsh.fltk.run()


mesh, ct, fm = gmshio.model_to_mesh(gmsh.model, comm, modelRank, gdim=3)
if mpiRank == modelRank:
    gmsh.finalize()

mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
elem = basix.ufl.element('Nedelec 1st kind H(curl)', mesh.basix_cell(), degree=2)
V = fem.functionspace(mesh, elem)
u = U.TrialFunction(V)
v = U.TestFunction(V)
normal = U.FacetNormal(mesh)
ds = U.Measure("ds", domain=mesh, subdomain_data = fm)
# Dirichlet BCs
facets = fm.find(1)
ubc = fem.Function(V)
ubc.interpolate(lambda x: np.array([0*x[0], 0*x[1], 0*x[2]], dtype=np.complex128))
dofs = fem.locate_dofs_topological(V, mesh.topology.dim-1, facets)
bc = fem.dirichletbc(ubc, dofs)
# Build rhs
uh = fem.Function(V, dtype=np.complex128)
kz = np.sqrt(k0 * k0 - (np.pi / a)**2)
f = ModesTE(kz, k0, a, b, 1, 0)
uh.interpolate(f.eval_H)
finc = 2.0j * k0 * eta0 * U.inner(U.cross(normal, uh), v) * ds(2)
uh.name = "H_inc"
Lb = 1.0j * kz * U.inner(U.cross(normal, u), U.cross(normal, v)) * ds(2)
Lc = 1.0j * kz * U.inner(U.cross(normal, u), U.cross(normal, v)) * ds(3)
t_start = time.perf_counter()

# Weak form
L = (U.inner(U.curl(u), U.curl(v)) - k0 * k0 * U.inner(u, v)) * U.dx
Df = fem.form(L+Lb+Lc)
Rf = fem.form(finc)

num_dofs_local = V.dofmap.index_map.size_local * V.dofmap.index_map_bs
num_dofs_global = V.dofmap.index_map.size_global * V.dofmap.index_map_bs
print(f"Number of dofs (owned) by rank {mpiRank}: {num_dofs_local}")
if mpiRank == modelRank:
    print(f"Number of dofs global: {num_dofs_global}")

def VecCreate(linear_form, bilin_form, BCs):
    vec = petsc.create_vector(linear_form)
    petsc.assemble_vector(vec, linear_form)
    petsc.apply_lifting(vec, [bilin_form], [[BCs]])
    vec.ghostUpdate(addv=P.InsertMode.ADD_VALUES, mode=P.ScatterMode.REVERSE)
    petsc.set_bc(vec, [BCs])
    vec.ghostUpdate(addv=P.InsertMode.INSERT, mode=P.ScatterMode.FORWARD)
    return vec

A = petsc.assemble_matrix(Df, bcs=[bc])
A.assemble()
b = VecCreate(Rf, Df, bc)

petsc_options = {
        "ksp_type": "preonly",
        "pc_type": "lu",
        "pc_factor_mat_solver": "mumps",
        "pc_factor_nonzeros_along_diagonal": 1e-8,
        "ksp_error_if_not_converged": False,
        "ksp_monitor": None }
opts = P.Options()
for key, value in petsc_options.items():
    opts[key] = value

solver = P.KSP().create(mesh.comm)
solver.setOperators(A)
solver.setFromOptions()
solver.setUp()
t_stop = time.perf_counter()

print("Solver setup time = {0}".format(t_stop - t_start))
sol = fem.Function(V)
sol.name="Solution"
solver.solve(b, sol.x.petsc_vec)
sol.x.petsc_vec.ghostUpdate(addv=P.InsertMode.INSERT, mode=P.ScatterMode.FORWARD)
t_solve = time.perf_counter()

print("Solve time = {0}".format(t_solve - t_stop))
Vw = basix.ufl.element('DG', mesh.basix_cell(), 0, shape=(mesh.geometry.dim, ))
W = fem.functionspace(mesh, Vw)
Et = fem.Function(W)
Et.interpolate(sol)
Et.name = "ElectricField"
with io.XDMFFile(mesh.comm, "Efield.xdmf", "w") as xx:
    xx.write_mesh(mesh)
    xx.write_function(Et)

sys.exit(0)

There could be a missing constant somewhere, but I would like to know if I am approaching the Sherman-Morrison method correctly using the PETSc vector manipulations. The results of the two methods should be equivalent.

2

I discovered a strange effect in the PETSc vector functions. Using the “PETSc.axpy(alpha, vec)” function seems to corrupt the solution vector. I tried

fsol.x.petsc_vec.zeroEntries()
fsol.x.petsc_vec.axpy(1.0, sol.x.petsc_vec)
print(fsol.x.petsc_vec.equal(sol.x.petsc_vec))

and it yields a “False” result when fsol should just be a duplicate of sol. What is happening here? Am I overlooking something? This is most likely corrupting my solution.