Good day eveyone. I am trying to work through a Spack installation of Fenicsx 0.9.0 with petsc 3.24.1 on a Ubuntu 24.04 system. The compilation process finishes without a hitch, but running an example crashes.
First the install parameters:
spack add petsc@3.24.1 -trilinos +complex +fortran +mumps +superlu-dist
spack add py-fenics-dolfinx@0.9.0 +petsc4py +slepc4py
spack add py-mpi4py@4.0.1
spack add py-numpy@2.1.2
spack add py-matplotlib
spack add py-pandas
spack add py-gmsh
spack install
I then run the following code (file with name Rect1.py) on one processor with the debugger
mpirun -np 1 python3 Rect1.py
and I see that the LU factorization stage fails immediately with a SEGV fault. Mesh building, matrix and vector assemblies execute without error. Is there some incompatibility between Dolfinx 0.9.0 and PETSc 3.24.1? Am I missing something in the compile configuration?
The code that I run:
# Example to test multi-mode waveguide absorbing boundary condition
# for a simple rectangular waveguide
import pdb
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.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)
def eval_E(self, x):
if((self.m != 0) and (self.n != 0)):
B = np.sqrt(self.n * self.m * 4.0 / (ZTE * self.a * self.b))
elif self.m == 0:
B = np.sqrt(self.n * 2.0 / (ZTE * self.a * self.b))
else:
B = np.sqrt(self.m * 2.0 / (ZTE * self.a * self.b))
ZTE = eta0 * self.k / self.beta
ey = -ZTE * B * np.sin((np.pi * self.m / self.a) * x[0]) * np.cos((np.pi * self.n / self.b) * x[1])
ex = ZTE * 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 = 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_type": "mumps",
"pc_factor_nonzeros_along_diagonal": 1e-8,
"ksp_error_if_not_converged": True,
"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()
pdb.set_trace()
solver.setUp()
print("***Made it to this point***")
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)
solver.view()
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)
The error message:
-> solver.setUp()
n
--------------------------------------------------------------------------
prterun noticed that process rank 0 with PID 66273 on node bandicoot exited on
signal 11 (Segmentation fault).
--------------------------------------------------------------------------
(Pdb)
Anyway, running the code modified for complex forms, I get an error on the entity map in the dolfinx.fem.form statement: