Thanks @Stein ! I took your advice and used MatCreateAIJ. The assembly worked well to build a nested matrix. However, the solve stage does not work. (Here is an example.)
# Example to test multi-mode waveguide absorbing boundary condition
# for a simple rectangular waveguide using nested PETSc matrices
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):
ZTE = self.k * eta0 / self.beta
if((self.m != 0) and (self.n != 0)):
B = np.sqrt(4.0 / (ZTE * self.a * self.b))
elif self.m == 0:
B = np.sqrt(2.0 / (ZTE * self.a * self.b))
else:
B = np.sqrt(2.0 / (ZTE * 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)
def eval_E(self, x):
ZTE = self.k * eta0 / self.beta
if((self.m != 0) and (self.n != 0)):
B = np.sqrt(4.0 * ZTE / (self.a * self.b))
elif self.m == 0:
B = np.sqrt(2.0 * ZTE / (self.a * self.b))
else:
B = np.sqrt(2.0 * ZTE / (self.a * self.b))
ey = -B * np.sin((np.pi * self.m / self.a) * x[0]) * np.cos((np.pi * self.n / self.b) * x[1])
ex = B * np.cos((np.pi * self.m / self.a) * x[0]) * np.sin((np.pi * self.n / self.b) * x[1])
ez = np.full_like(hx, 0.0+0j, dtype=np.complex128)
return(ex, ey, ez)
eps = 1.0e-4
a = 22.86 #WG width
b = 10.16 #WG height
l = 25.8 # WG length
f = 10.5 #frequency of operation GHz
c = 300.0 #speed of light
lm = 8.20
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', 8.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)
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
P1 = mesh.comm.allreduce(fem.assemble_scalar(fem.form(U.inner(ua, ua) * ds(2))) * Zm , op=M.SUM)
P2 = mesh.comm.allreduce(fem.assemble_scalar(fem.form(U.inner(ua, ua) * ds(3))) * Zm , op=M.SUM)
print(P1, P2)
t_start = time.perf_counter()
# Weak form
L = (U.inner(U.curl(u), U.curl(v)) - k0 * k0 * U.inner(u, v)) * U.dx
Rf = fem.form(finc)
Df = fem.form(L)
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
def VecCreate(linear_form, bilin_form, BCs):
vec = petsc.create_vector(linear_form)
with vec.localForm() as lb:
lb.set(0.0)
vec.set(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.assemble_matrix(Df, bcs=[bc])
A.assemble()
b = VecCreate(Rf, Df, bc)
u1 = VecCreate(um1a, Df, bc)
u2 = VecCreate(um1b, Df, bc)
sol = fem.Function(V)
U = P.Mat().createAIJ(size=(b.size, 2))
U.setUp()
print("Vector size ={0}".format(b.getSize()), flush=True)
print("MPI process ownership ranges" + str(b.getOwnershipRanges()), flush=True)
iz = b.getOwnershipRanges()
nnb = 0
for ii, index in enumerate(iz):
if ii != 0:
if mpiRank == ii-1:
for jj in range(iz[ii-1], index):
aa = b.getValue(jj)
if(np.abs(aa) > 0.01*eps):
print("Proc ={0}, index={1}, Entry={2}".format(ii, jj, aa), flush=True)
nnb += 1
print(mpiRank, nnb, flush=True)
print("Vector size ={0}".format(u1.getSize()), flush=True)
print("MPI process ownership ranges" + str(u1.getOwnershipRanges()), flush=True)
iz = u1.getOwnershipRanges()
nnu1 = 0
for ii, index in enumerate(iz):
if ii != 0:
if mpiRank == ii-1:
for jj in range(iz[ii-1], index):
aa = u1.getValue(jj)
if(np.abs(aa) > 0.001*eps):
U.setValue(jj, 0, 1j * kz * aa, addv=P.InsertMode.INSERT_VALUES)
print("Proc ={0}, index={1}, Entry={2}".format(ii, jj, 1j * kz * aa), flush=True)
nnu1 += 1
print(mpiRank, nnu1, flush=True)
print("Vector size ={0}".format(u2.getSize()), flush=True)
print("MPI process ownership ranges" + str(u2.getOwnershipRanges()), flush=True)
iz = u2.getOwnershipRanges()
nnu2 = 0
for ii, index in enumerate(iz):
if ii != 0:
if mpiRank == ii-1:
for jj in range(iz[ii-1], index):
aa = u2.getValue(jj)
if(np.abs(aa) > 0.001*eps):
U.setValue(jj, 1, 1j * kz * aa, addv=P.InsertMode.INSERT_VALUES)
print("Proc ={0}, index={1}, Entry={2}".format(ii, jj, 1j * kz * aa), flush=True)
nnu2 += 1
print(mpiRank, nnu2, flush=True)
U.assemble()
Ut = U.transpose()
# Identity block
Id = P.Mat().createAIJ(size=(2, 2))
Id.setUp()
Id.setValue(0, 0, 1.0, addv=P.InsertMode.INSERT_VALUES)
Id.setValue(1, 1, 1.0, addv=P.InsertMode.INSERT_VALUES)
Id.assemble()
A_new = P.Mat()
A_new.createNest([[A, U], [Ut, Id]], comm = comm)
A_new.setUp()
A_new.assemble()
b_0 = P.Vec().create(comm=comm)
b_0.setSizes(2)
b_0.setUp()
b_0.assemble()
print("***********Last line!*************")
x_0 = P.Vec().create(comm=comm)
x_0.setSizes(2)
x_0.setUp()
x_0.assemble()
b_new = P.Vec()
b_new.createNest([b, b_0], comm = comm)
b_new.setUp()
b_new.assemble()
x_new = P.Vec()
x_new.createNest([sol.x.petsc_vec, x_0], comm=comm)
x_new.setUp()
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_new)
solver.setFromOptions()
solver.setUp()
t_stop = time.perf_counter()
solver.solve(b_new, x_new)
print(A_new.sizes)
sys.exit(0)
I get this error:
petsc4py.PETSc.Error: error code 92
[2] KSPSetUp() line 406 in ./src/ksp/ksp/interface/itfunc.c
[2] PCSetUp() line 1017 in ./src/ksp/pc/interface/precon.c
[2] PCSetUp_LU() line 82 in ./src/ksp/pc/impls/factor/lu/lu.c
[2] MatGetFactor() line 4756 in ./src/mat/interface/matrix.c
[2] See https://petsc.org/release/overview/linear_solve_table/ for possible LU and Cholesky solvers
[2] Could not locate a solver type for factorization type LU and matrix type nest.
solver.setUp()
I have seen from other posts that MUMPS LU solver does not work for nested matrices.
This leads me to a couple of questions.
- Is there a way to convert the nested form into a form suitable for LU factorization solvers?
- Would it be easier to avoid the nested construction? What type of matrix is suitable for direct methods? How would I go about creating a monolithic sparse matrix for use with an LU solver?
Unfortunately, I cannot find much documentation on the topic. Most assumes that iterative solvers are used (Navier Stokes p-v formulations, etc.) I need to manually insert rows/columns in the PETSc matrix that do not correspond to FEM degrees of freedom, but represent external constraints.
. Since I am not a particularly good computer administrator, what is the best way to force the installation of a newer version of petsc4py (3.20 or newer) ? Even pip wants to install v. 3.19.