I am trying to solve a generalised eigenvalue problem Aq=\lambda Bq with boundary conditions on q.
Here’s a toy problem, with a 2D Laplacian and homogeneous Dirichlet conditions on the unit square. Theoretical solution for the eigenvalues should be -\pi^2(n^2+m^2) but a nave approach spits them wrong :
import numpy as np
from slepc4py import SLEPc as slp
from mpi4py.MPI import COMM_WORLD
mesh=dolfinx.UnitSquareMesh(COMM_WORLD,4,4)
FE=ufl.FiniteElement("Lagrange",mesh.ufl_cell(),1)
V=dolfinx.FunctionSpace(mesh,FE)
u=ufl.TrialFunction(V)
v=ufl.TestFunction(V)
# Defining the forms
a=ufl.inner(ufl.grad(u),ufl.grad(v))*ufl.dx
b=ufl.inner( u, v )*ufl.dx
# Boundary conditions
zeros=dolfinx.Function(V)
with zeros.vector.localForm() as zero_loc: zero_loc.set(0)
dim = mesh.topology.dim - 1
# Create facet to cell connectivity required to determine boundary facets
mesh.topology.create_connectivity(dim, mesh.topology.dim)
boundary = np.where(np.array(dolfinx.cpp.mesh.compute_boundary_facets(mesh.topology)) == 1)[0]
boundary_dofs = dolfinx.fem.locate_dofs_topological(V, dim, boundary)
BC = dolfinx.DirichletBC(zeros, boundary_dofs)
# Assembling matrices
A=dolfinx.fem.assemble_matrix(a, bcs=[BC]); A.assemble()
B=dolfinx.fem.assemble_matrix(b, bcs=[BC]); B.assemble() # Presumed location of problem
# Solver
E = slp.EPS(); E.create(COMM_WORLD)
E.setOperators(A,B) # Solve A*x=lambda*B*x
E.setProblemType(slp.EPS.ProblemType.GHEP)
E.setWhichEigenpairs(E.Which.SMALLEST_MAGNITUDE)
E.setDimensions(10)
E.solve()
for i in range(10):
print(np.abs(np.pi**2*((i//2)**2+(i//2+i%2)**2)+E.getEigenvalue(i)))
Output on my end (Docker installation, dolfinx in complex mode) :
Quick update. I’ve moved to scipy for easier matrix slicing and tried my luck at axing away the homogeneous dofs using this post to obtain total DoFs and taking its complementary to obtain a mask of free indices and that post to transform the PETSc Mat into a scipy CSR :
import dolfinx, ufl
import numpy as np
import scipy.sparse as sps
import scipy.sparse.linalg as la
from mpi4py.MPI import COMM_WORLD
n=100
mesh=dolfinx.UnitSquareMesh(COMM_WORLD,n,n)
FE=ufl.FiniteElement("Lagrange",mesh.ufl_cell(),1)
V=dolfinx.FunctionSpace(mesh,FE)
u=ufl.TrialFunction(V)
v=ufl.TestFunction(V)
# Defining the forms
a=ufl.inner(ufl.grad(u),ufl.grad(v))*ufl.dx
b=ufl.inner( u, v )*ufl.dx
# Identify BC DoFs
dim = mesh.topology.dim - 1
mesh.topology.create_connectivity(dim, mesh.topology.dim)
boundary = np.where(np.array(dolfinx.cpp.mesh.compute_boundary_facets(mesh.topology)) == 1)[0]
boundary_dofs = dolfinx.fem.locate_dofs_topological(V, dim, boundary)
# Total number of DoFs
N = V.dofmap.index_map.size_global * V.dofmap.index_map_bs
# Free DoFs
freeinds = np.setdiff1d(range(N),boundary_dofs,assume_unique=True).astype(np.int32)
# Matrices
A = dolfinx.fem.assemble_matrix(a); A.assemble()
# Convert from PeTSc to Scipy for easy slicing
ai, aj, av = A.getValuesCSR()
# Slice away BCs
A = sps.csr_matrix((av, aj, ai))[freeinds,:][:,freeinds]
# Matrices
B = dolfinx.fem.assemble_matrix(b); B.assemble()
bi, bj, bv = B.getValuesCSR()
B = sps.csr_matrix((bv, bj, bi))[freeinds,:][:,freeinds]
# Solver
k=5
vals, vecs = la.eigs(A, k=k, M=B, sigma=0)
for i in range(k): print(np.abs(np.pi**2*((i//2)**2+(i//2+i%2)**2)+vals[i]))
I still get absurd errors, starting at 20 and going all the way up to 177. Anyone has an idea out there ?
I’ve tried seomthing even more brutal with imposing some rows of A to be ones and zeros with corresponding rows oof B being 0 to no avail
# Total number of DoFs
N = V.dofmap.index_map.size_global * V.dofmap.index_map_bs
n_bc=boundary_dofs.size
BC_id,BC_zeros=np.zeros((n_bc,N)),np.zeros((n_bc,N))
for i in range(n_bc):
BC_id[i,boundary_dofs[i]]=1
# Matrices
A = dolfinx.fem.assemble_matrix(a); A.assemble()
# Convert from PeTSc to Scipy for easy slicing
ai, aj, av = A.getValuesCSR()
# Slice away BCs
A = sps.csr_matrix((av, aj, ai))
A[boundary_dofs,:]=BC_id
# Same for B
B = dolfinx.fem.assemble_matrix(b); B.assemble()
bi, bj, bv = B.getValuesCSR()
B = sps.csr_matrix((bv, bj, bi))
B[boundary_dofs,:]=BC_zeros
Ok sorry for the spam I had my theoretical values wrong. The following code works, and it’s the most elegant / efficient way I’ve found to impose homogeneous BCs for a generalised eigenvalue pb.
Any suggestions about how to do the same thing in full PETSc or efficiently for an arbitrary BC is welcome.
import dolfinx, ufl
import numpy as np
import scipy.sparse as sps
import scipy.sparse.linalg as la
from mpi4py.MPI import COMM_WORLD
n=100
mesh=dolfinx.UnitSquareMesh(COMM_WORLD,n,n)
FE=ufl.FiniteElement("Lagrange",mesh.ufl_cell(),1)
V=dolfinx.FunctionSpace(mesh,FE)
u=ufl.TrialFunction(V)
v=ufl.TestFunction(V)
# Defining the forms
a=ufl.inner(ufl.grad(u),ufl.grad(v))*ufl.dx
b=ufl.inner( u, v )*ufl.dx
# Identify BC DoFs
dim = mesh.topology.dim - 1
mesh.topology.create_connectivity(dim, mesh.topology.dim)
boundary = np.where(np.array(dolfinx.cpp.mesh.compute_boundary_facets(mesh.topology)) == 1)[0]
boundary_dofs = dolfinx.fem.locate_dofs_topological(V, dim, boundary)
# Total number of DoFs
N = V.dofmap.index_map.size_global * V.dofmap.index_map_bs
# Free DoFs
freeinds = np.setdiff1d(range(N),boundary_dofs,assume_unique=True).astype(np.int32)
# Matrices
A = dolfinx.fem.assemble_matrix(a); A.assemble()
# Convert from PeTSc to Scipy for easy slicing
ai, aj, av = A.getValuesCSR()
# Slice away BCs
A = sps.csr_matrix((av, aj, ai))[freeinds,:][:,freeinds]
# Same for B
B = dolfinx.fem.assemble_matrix(b); B.assemble()
bi, bj, bv = B.getValuesCSR()
B = sps.csr_matrix((bv, bj, bi))[freeinds,:][:,freeinds]
# Solver
k=5
vals, vecs = la.eigs(A, k=k, M=B, sigma=0)
ths=np.pi**2*np.array([2,5,5,8,10])
print("Theoretical value | Calculated one | Difference")
for i in range(k): print(f"{ths[i]:1.14f}"+" | "+f"{np.real(vals[i]):.4f}"+f"{np.imag(vals[i]):+.0e}"+"j | "+f"{np.abs(ths[i]-vals[i]):00.2e}")
As a follow up question, linked to this post, why is locate_dofs_geometrical returning a list of two arrays ? I expected a single 1D arrray of integers pointing to dofs in a direct manner, going from 0 to self.Space.dofmap.index_map.size_global * self.Space.dofmap.index_map_bs-1 and identifying each dof in a unique manner.
Now I’m unsure my above hack works as expected, given freeinds has been computed naively from a set difference.