Dea All,
I tried to solve a general eigenvalue problem and got the error frequently
petsc4py.PETSc.Error: error code 71
[0] EPSSolve() line 136 in /usr/local/slepc/src/eps/interface/epssolve.c
[0] EPSSetUp() line 350 in /usr/local/slepc/src/eps/interface/epssetup.c
[0] STSetUp() line 582 in /usr/local/slepc/src/sys/classes/st/interface/stsolve.c
[0] STSetUp_Shift() line 107 in /usr/local/slepc/src/sys/classes/st/impls/shift/shift.c
[0] KSPSetUp() line 406 in /usr/local/petsc/src/ksp/ksp/interface/itfunc.c
[0] PCSetUp() line 1015 in /usr/local/petsc/src/ksp/pc/interface/precon.c
[0] PCSetUp_LU() line 131 in /usr/local/petsc/src/ksp/pc/impls/factor/lu/lu.c
[0] MatLUFactorNumeric() line 3195 in /usr/local/petsc/src/mat/interface/matrix.c
[0] MatLUFactorNumeric_SeqAIJ_Inode() line 1725 in /usr/local/petsc/src/mat/impls/aij/seq/inode.c
[0] MatPivotCheck() line 829 in /usr/local/petsc/include/petsc/private/matimpl.h
[0] MatPivotCheck_none() line 812 in /usr/local/petsc/include/petsc/private/matimpl.h
[0] Zero pivot in LU factorization: https://www.mcs.anl.gov/petsc/documentation/faq.html#zeropivot
[0] Zero pivot row 3422 value 0. tolerance 2.22045e-14
Is there any method to solve this problem? My code is
import pygmsh
import numpy as np
import gmsh, meshio
from mpi4py import MPI
from slepc4py import SLEPc
from ufl import curl, dx, FiniteElement, grad, inner, Measure, TestFunctions, TrialFunctions
from dolfinx import Constant, Function, FunctionSpace, geometry, cpp, DirichletBC
from dolfinx.fem import locate_dofs_topological, assemble_matrix
from scipy.constants import c, epsilon_0
from dolfinx.io import XDMFFile
from petsc4py import PETSc
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
if prune_z:
out_mesh.prune_z_0()
return out_mesh
with pygmsh.occ.Geometry() as geom:
core_r = 0.0015
max_r = 4000
inner_lcar = 0.0001
outer_lcar = 50
theta = np.pi / 8
p0 = [0, 0]
p1 = [core_r, 0]
p2 = [core_r * np.cos(theta), core_r * np.sin(theta)]
p3 = [max_r * np.cos(theta), max_r * np.sin(theta)]
p4 = [max_r, 0]
point_coords = [p0, p1, p2, p3, p4]
lcars = [inner_lcar, inner_lcar, inner_lcar, outer_lcar, outer_lcar]
points = []
for p, lcar in zip(point_coords, lcars):
points.append(geom.add_point([p[0], p[1]], lcar))
line1 = geom.add_line(points[0], points[1])
inner_arc = geom.add_circle_arc(points[1], points[0], points[2])
line2 = geom.add_line(points[2], points[0])
inner_loop_curves = [line1, inner_arc, line2]
line3 = geom.add_line(points[2], points[3])
outer_arc = geom.add_circle_arc(points[3], points[0], points[4])
line4 = geom.add_line(points[4], points[1])
outer_loop_curves = [inner_arc, line3, outer_arc, line4]
inner_curve_loop = geom.add_curve_loop(inner_loop_curves)
outer_curve_loop = geom.add_curve_loop(outer_loop_curves)
inner_surface = geom.add_plane_surface(inner_curve_loop)
outer_surface = geom.add_plane_surface(outer_curve_loop)
geom.add_physical(outer_arc, 'Dirichlet')
geom.add_physical(inner_surface, 'Conductor')
geom.add_physical(outer_surface, 'Air')
mesh = geom.generate_mesh(2)
# Read in mesh
gmsh.write("mesh.msh")
msh = meshio.read("mesh.msh")
# Create and save one file for the mesh, and one file for the facets
triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
line_mesh = create_mesh(msh, "line", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mt.xdmf", line_mesh)
with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf_infile:
mesh = xdmf_infile.read_mesh(name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim)
mvc_subdomain = xdmf_infile.read_meshtags(mesh, "Grid")
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
with XDMFFile(MPI.COMM_WORLD, "mt.xdmf", "r") as xdmf_infile:
mvc_boundaries = xdmf_infile.read_meshtags(mesh, "Grid")
V = FiniteElement('N1curl', mesh.ufl_cell(), 2)
Q = FiniteElement('P', mesh.ufl_cell(), 3)
W = FunctionSpace(mesh, V*Q)
u, p = TrialFunctions(W)
v, q = TestFunctions(W)
omega = 1e4
sigma_core = 8678
sigma_corona = 5e-6
k0_sqr = (omega/c)**2
S = FunctionSpace(mesh, ('DG', 0))
e_r = Function(S)
x = S.tabulate_dof_coordinates()
for i in range(x.shape[0]):
if mvc_subdomain.values[i] == 2:
e_r.vector.setValueLocal(i, 1-1j*sigma_core/(omega*epsilon_0))
else:
e_r.vector.setValueLocal(i, 1)
theta_sqr = k0_sqr * 1.1
print(theta_sqr)
a = 1/theta_sqr*(inner(curl(u), curl(v)) - k0_sqr*e_r*inner(u, v))*dx
b = inner(grad(p), v)*dx + inner(u, v)*dx + inner(grad(p), grad(q))*dx + inner(u, grad(q))*dx - k0_sqr*e_r*inner(p, q)*dx
c = b + a
electric_wall = Function(W)
with electric_wall.vector.localForm() as bc_local:
bc_local.set(0)
bndry_facets = np.where(mvc_boundaries.values == 0)[0]
bdofs = locate_dofs_topological(W, 1, bndry_facets)
bc = DirichletBC(electric_wall, bdofs)
B = assemble_matrix(b, [bc])
B.assemble()
C = assemble_matrix(c, [bc])
C.assemble()
eps = SLEPc.EPS()
eps.setOperators(B, C)
eps.setDimensions(50)
eps.solve()
n = eps.getConverged()
print(n)
xr, xi = B.getVecs()
for j in range(n):
lam = eps.getEigenpair(j, xr, xi)
if (((1-1/lam)*theta_sqr)**0.5).real > 0.0005 and abs((((1-1/lam)*theta_sqr)**0.5).imag) < 1:
print(lam)
print(((1-1/lam)*theta_sqr)**0.5, 'xxx')
x_r = Function(W)
xr.copy(x_r.vector)
x_i = Function(W)
xi.copy(x_i.vector)
xt_r, xz_r = x_r.split()
xt_i, xz_i = x_i.split()
with XDMFFile(MPI.COMM_WORLD, "Et_r.xdmf", 'w',
encoding=XDMFFile.Encoding.HDF5) as file:
file.write_mesh(mesh)
file.write_function(xt_r)
with XDMFFile(MPI.COMM_WORLD, "Ez_r.xdmf", 'w',
encoding=XDMFFile.Encoding.HDF5) as file:
file.write_mesh(mesh)
file.write_function(xz_r)
with XDMFFile(MPI.COMM_WORLD, "Et_i.xdmf", 'w',
encoding=XDMFFile.Encoding.HDF5) as file:
file.write_mesh(mesh)
file.write_function(xt_i)
with XDMFFile(MPI.COMM_WORLD, "Ez_i.xdmf", 'w',
encoding=XDMFFile.Encoding.HDF5) as file:
file.write_mesh(mesh)
file.write_function(xz_i)
I’m grateful for your help.