Dear Fenics community, inspired by the post Does Fenicsx support spectral element methods? , I built an attempt to apply spectral element methods to the computation of an eigenvalue problem (specifically, the quasi-normal modes of a Kerr black-hole).
The code below works perfectly. By using a one-mesh cell and setting a high enough polynomial degree (> 10), it returns (with the set of parameters I choose here) something very close to the (theoretical result for the) fundamental quasi-normal mode omega = 0.74734 - 0.17793*1j .
As stated in the link above, one cell + high polynomial degree is equivalent to applying a spectral method, while increasing the number of cells should be equivalent to applying a spectral element method.
My issue is this problem seems to be solvable by spectral methods only, not spectral element methods. Whenever I increase the number of cells to anything beyond one (even without decreasing my polynomial degree), the accuracy plummets. For example, with 1 cell and degree 10, the results are accurate up to the 10^-5 decimal place. But with the same degree (or even higher) and 2 cells, the results come out with a relative error of > 10%. It looks like the problem fundamentally changes when I increase the number of cells, because even if I keep increasing the polynomial degree, I do not recover accuracy: the eigenvalues obtained converge to a given value with the error stated above.
I suspect there is an issue in the boundary between elements, since the eigenfunctions recovered look discontinuous there, but I am no expert on the topic. Therefore, I would be extremely grateful if anyone could point out why it is not possible to solve an eigenvalue problem like this via spectral element methods or, in case it is, pinpoint what I might be doing wrong.
Thanks in advance 
import numpy as np
from math import pi
import ufl
from mpi4py import MPI
from dolfinx.fem.petsc import assemble_matrix
from dolfinx import mesh
from dolfinx.fem import functionspace
from dolfinx import fem
import basix
from basix import CellType, ElementFamily, LagrangeVariant, QuadratureType
from basix.ufl import custom_element
from dolfinx.fem.petsc import assemble_matrix
from slepc4py import SLEPc
import numpy as np
from petsc4py import PETSc
import time
######### PHYSICAL PARAMETERS #########
M_mass = 0.5
alpha = 0 # a/M, in (0,1)
a_rot = alpha*M_mass # a, in (0,M)
r_p = M_mass*(1 + np.sqrt(1 - (a_rot**2)/(M_mass**2) ) ) # Outer horizon
kappa = alpha/(1 + np.sqrt(1 - alpha**2)) # Dimensionless angular momentum parameter
#### Gauge choice ####
lambda_param = r_p
######## Spin and azimuthal mode #######
p_mode = -2
m_mode = 1
####### Angular equation eigenvalue #######
A_lm = 4
######### SPECTRAL ELEMENT AND DOMAIN DEFINITION #########
N_x = 1
poly_degree = 10
# 1D Problem === Interval type cell
cell_type = CellType.interval
# Element creation using Gauss-Lobatto-Legendre (GLL) nodes
nodal_element = basix.create_element(ElementFamily.P,
cell_type,
poly_degree,
LagrangeVariant.gll_isaac)
ufl_nodal_element = custom_element(cell_type,
[],
nodal_element.wcoeffs, # pyright: ignore
nodal_element.x, # pyright: ignore
nodal_element.M, # pyright: ignore
nodal_element.interpolation_nderivs,
nodal_element.map_type,
nodal_element.sobolev_space,
nodal_element.discontinuous,
nodal_element.degree,
nodal_element.embedded_subdegree,
nodal_element.polyset_type)
t_0 = time.time()
# Meshing
msh = mesh.create_interval(MPI.COMM_WORLD, N_x, [0.0, 1.0])
tdim = msh.topology.dim
fdim = tdim - 1
msh.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(msh.topology)
# Function spaces
V_space = functionspace(msh, ufl_nodal_element)
V_space_plot = functionspace(msh, ("Lagrange", poly_degree))
# Defining the trial and test functions
u_trial = ufl.TrialFunction(V_space)
v_main = ufl.TestFunction(V_space)
######### CUSTOM QUADRATURE RULE #########
pts, wts = basix.make_quadrature(cell_type,
2*poly_degree-1,
rule=QuadratureType(2))
metadata = {"quadrature_rule": "custom",
"quadrature_points": pts, "quadrature_weights": wts}
dx_meas = ufl.dx(domain=msh, metadata=metadata)
######### TEUKOLSKY COEFFICIENT FUNCTIONS #########
def m_0(x):
term_kerr = - kappa**2 - 4*(kappa**2)*(x**2)*(kappa**2 + 1) - 4*(kappa**4)*(x**2)*(kappa**2 + 1)
term_kerr += + 4*(kappa**4)*x*(kappa**2 + 1)
return -(4*x*(kappa**2 + 1)**2 + 4*(kappa**2 + 1)**2 + term_kerr )
def c_1(x):
kerr_term = -6*(kappa**2)*(x**2) + 4*(kappa**4)*(x**2)*(x-1)
return (kappa**2*x -1)*(4*(x**2)) + 2 + kerr_term
def c_0(x):
part_1 = 2*x*(kappa**2 + 1) - 2*kappa**2*(x**2) - 2*(kappa**4)*(x**2) + 2*kappa*m_mode*1j
part_2 = 2*x*(kappa**2 + 1)**2 + 2*(kappa**4)*x - 2*p_mode*(kappa**2 + 1) - 4*(kappa**2)*(x**2)*(kappa**2 + 1)
part_3 = 2*x*(kappa**2 + 1)*(p_mode*(kappa**2 + 1) + 2*kappa*m_mode*1j)
return -part_1-part_2-part_3
def k_2(x):
return (kappa**2*(x) - 1)*(x**2)*(x-1)
def k_1(x):
sig1 = p_mode*(kappa**2 + 1)
part_1 = 2*x + 2*p_mode*x
part_2 = -(x**2)*(3*kappa**2 - 4*kappa**2*x + sig1 + 3 + 2*kappa*m_mode*1j)
return part_1 + part_2
def k_0(x):
sig1 = p_mode*(kappa**2 + 1)
return -(A_lm + x*(kappa**2 + 1) -2*(kappa**2)*(x**2) + x*(sig1 + 2*kappa*m_mode*1j))
# Coefficient interpolation into the function space
m_0_term = fem.Function(V_space)
m_0_term.interpolate(lambda x: m_0(x[0]))
c_1_term = fem.Function(V_space)
c_1_term.interpolate(lambda x: c_1(x[0]))
c_0_term = fem.Function(V_space)
c_0_term.interpolate(lambda x: c_0(x[0]))
k_2_term = fem.Function(V_space)
k_2_term.interpolate(lambda x: k_2(x[0]))
k_1_term = fem.Function(V_space)
k_1_term.interpolate(lambda x: k_1(x[0]))
k_0_term = fem.Function(V_space)
k_0_term.interpolate(lambda x: k_0(x[0]))
dev_k2 = fem.Function(V_space)
dev_k2.interpolate(lambda x: -x[0]*(-4*(kappa**2)*(x[0]**2) + 3*(kappa**2)*x[0] + 3*x[0] - 2))
######### WEAK FORMULATION AND MATRIX ASSEMBLY #########
#### M Matrix
M_eq = m_0_term*ufl.inner(u_trial,v_main)*dx_meas
#### C Matrix
C_eq = c_1_term*ufl.inner(u_trial.dx(0),v_main)*dx_meas + c_0_term*ufl.inner(u_trial,v_main)*dx_meas
#### K Matrix
K_eq = - k_2_term*ufl.inner(u_trial.dx(0),v_main.dx(0))*dx_meas
K_eq += - dev_k2*ufl.inner(u_trial.dx(0),v_main)*dx_meas + k_1_term*ufl.inner(u_trial.dx(0),v_main)*dx_meas
K_eq += k_0_term*ufl.inner(u_trial,v_main)*dx_meas
M_mat = assemble_matrix(fem.form(M_eq))
M_mat.assemble()
C_mat = assemble_matrix(fem.form(C_eq))
C_mat.assemble()
K_mat = assemble_matrix(fem.form(K_eq))
K_mat.assemble()
A = [ ]
A.append(K_mat)
A.append(C_mat)
A.append(M_mat)
print("Problem assembled successfully... now solving")
######### SLEPc SOLVER CONFIGURATION (PEP) #########
target = -1j*lambda_param*(1)
nev = 4
tol = 1e-7
max_it = 1000
comm = MPI.COMM_WORLD
pep = SLEPc.PEP().create(comm)
pep.setOperators(A)
pep.setProblemType(SLEPc.PEP.ProblemType.GENERAL)
pep.setDimensions(nev=nev, ncv=max(15, 2*nev + 10), mpd=max(15, 2*nev + 10))
pep.setType(SLEPc.PEP.Type.TOAR)
st = pep.getST()
st.setType(SLEPc.ST.Type.SINVERT)
ksp = st.getKSP()
ksp.setType(PETSc.KSP.Type.PREONLY)
pep.setRefine(ref = SLEPc.PEP.Refine.SIMPLE, npart = 1, tol = 0.1*tol, its= 3, scheme= SLEPc.PEP.RefineScheme.SCHUR)
print("SCHUR Refine activated")
pep.setTarget(target)
pep.setWhichEigenpairs(SLEPc.PEP.Which.TARGET_MAGNITUDE)
pep.setTolerances(tol, max_it)
# Solve problem
pep.solve()
elapsed = time.time() - t_0
print(f"Calculation completed. Problem solved in {elapsed:.4f} s")
nconv = pep.getConverged()
if comm.rank == 0:
print(f"Converged eigenvalues: {nconv} QNM candidates.")
physical_eigvals = []
efuns = fem.Function(V_space)
efuns_list = []
error_list = []
for i in range(nconv):
eigval = pep.getEigenpair(i, efuns.x.petsc_vec)
error = pep.computeError(i)
if abs(eigval.real) < 10:
error_list.append(error)
efun_vals = efuns.x.array
efuns_list.append(efun_vals)
physical_eigvals.append(eigval/(-1j*lambda_param))
######### EXTRACTION AND POST-PROCESSING OF RESULTS #########
if MPI.COMM_WORLD.rank == 0:
print("\n--- Sorted Quasinormal Modes (QNMs) ---")
print(f"Found {len(physical_eigvals)} physical eigenpairs after sorting.")
for i, eigval in enumerate(physical_eigvals):
print(f"QNM {i}: {2*M_mass*eigval.real:.6e} + {2*M_mass*eigval.imag:.6e}i, error {error_list[i]}")