Dear Fenics users, I am trying to obtain the spectrum of eigenvalues ω for Regge-Wheeler’s equation (Wikipedia: Regge, Wheeler and Zerilli equations) with Neumann boundary conditions:
where
but the task is proving to be too challenging to tackle using SLEPc (as in Electromagnetic modal analysis for a waveguide).
In order to try a different approach, I tried to input one of the expected eigenvalues of the PDE and simply solve it using dolfinx. This works perfectly (code attached at the end of the post). But it requires me to know the eigenvalue in advance. Is there a way to work it out inside-out i.e., retrieving the value of ω knowing the PDE, boundary conditions and some sort of normalization condition? I believe it should be possible, but the fact that the BCs depend on said parameter is proving this task difficult (main reason why SLEPc does not seem to work)
Thank you very much in advance. This is my code:
from dolfinx import default_scalar_type
from dolfinx.fem import (Constant, Function, functionspace, assemble_scalar,
dirichletbc, form, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_unit_square, locate_entities, meshtags
from dolfinx.plot import vtk_mesh
from dolfinx import mesh
import ufl
from petsc4py import PETSc
#Code must be run with DOLFINx (PETSc) complex mode.
from mpi4py import MPI
from ufl import (FacetNormal, Measure, SpatialCoordinate, TestFunction, TrialFunction,
div, dot, dx, grad, inner, lhs, rhs)
import scipy
import numpy as np
import matplotlib.pyplot as plt
if not np.issubdtype(default_scalar_type, np.complexfloating):
print("Code should only be executed with DOLFINx complex mode")
exit(0)
M = 0.5
ell = 2
l = ell
w0 = 0.5*(0.747343 - 0.177925j)/M
s0 = 1j*w0
def rstar(r):
return r + 2*M*np.log(r/(2*M) - 1)
r_min = rstar(2.002*M)
r_max = rstar(80*M)
n_elements = 1000
domain = mesh.create_interval(MPI.COMM_WORLD, n_elements, [r_min, r_max])
mesh = domain
x = SpatialCoordinate(mesh)
# Define physical parameters and boundary condtions
n = FacetNormal(mesh)
kappa = Constant(mesh, default_scalar_type(1))
rplus = Constant(mesh, default_scalar_type(np.complex128(s0)))
rminus = Constant(mesh, default_scalar_type(-np.complex128(s0)))
s = Constant(mesh, default_scalar_type(0))
# Define function space and standard part of variational form
V = functionspace(mesh, ("Lagrange", 1))
u, v = TrialFunction(V), TestFunction(V)
def omegaux(x):
return scipy.special.wrightomega(-1 + x/(2*M))
def r(x):
return 2*M*(omegaux(x) + 1)
def f(x):
return 1 - 2*M/r(x)
V_l = Function(V)
V_l.interpolate(lambda x: f(x[0]) * (l * (l + 1) / (r(x[0])**2) - 6 * M / (r(x[0])**3)) )
F = kappa*inner(grad(u), grad(v))*dx + inner(V_l*u,v)*dx + (s0**2)*inner(u,v)*dx
u_a = lambda x: ufl.exp(x[0]*w0)
u_b = lambda x: ufl.exp(-x[0]*w0)
x = SpatialCoordinate(mesh)
# Define physical parameters and boundary condtions
boundaries = [(1, lambda x: np.isclose(x[0], r_min)),
(2, lambda x: np.isclose(x[0], r_max))]
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)
class BoundaryCondition():
def __init__(self, type, marker, values):
self._type = type
if type == "Dirichlet":
u_D = Function(V)
u_D.interpolate(values)
facets = facet_tag.find(marker)
dofs = locate_dofs_topological(V, fdim, facets)
self._bc = dirichletbc(u_D, dofs)
elif type == "Neumann":
self._bc = inner(values, v) * ds(marker)
elif type == "Robin":
self._bc = values[0] * inner(u-values[1], v)* ds(marker)
else:
raise TypeError("Unknown boundary condition: {0:s}".format(type))
@property
def bc(self):
return self._bc
@property
def type(self):
return self._type
# Define the condition
n = FacetNormal(mesh)
g1 = -dot(n, grad(u_a(x)))
g2 = -dot(n, grad(u_b(x)))
A = 1+1j
ufun = Function(V)
ufun.interpolate(lambda x: x[0]-x[0]+A)
u_bc = Function(V)
u_bc.interpolate(ufun)
boundary_conditions = [BoundaryCondition("Neumann", 1, g2),
BoundaryCondition("Neumann", 2, g1)]
bcs = []
for condition in boundary_conditions:
if condition.type == "Dirichlet":
bcs.append(condition.bc)
else:
F += condition.bc
# Solve linear variational problem
a = lhs(F)
L = rhs(F)
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
# Visualize solution
V_plot = (M**2)*V_l.x.array
rstar_plot = np.linspace(r_min, r_max, n_elements+1)
def r_from_tortoise(r_star):
# Lambert W function implementation for r(r*)
w = scipy.special.lambertw(np.exp((r_star - r_min)/(2*M) - 1))
return 2*M*(1 + w.real) # Taking real part
r_plot = r_from_tortoise(rstar_plot)
u_plotreal = uh.x.array.real
u_plotimag = uh.x.array.imag
max_amp = max(np.max(np.abs(u_plotreal)), np.max(np.abs(u_plotimag)))
plt.figure(figsize=(14, 7))
plt.rcParams.update({'font.size': 14})
plt.subplot(1, 2, 1)
plt.plot(r_plot, V_plot)
plt.title('Potencial de Regge-Wheeler V(r*)', fontsize=16)
plt.xlabel('r*', fontsize=15)
plt.ylabel(r'$V(r^*) \cdot M^2$', fontsize=15) # Usamos notación LaTeX para V(r*) * M^2
plt.grid(True)
plt.tick_params(axis='both', which='major', labelsize=13)
plt.subplot(1, 2, 2)
plt.plot(r_plot, u_plotreal/max_amp, label='Re(Ψ)')
plt.plot(r_plot, u_plotimag/max_amp, label='Im(Ψ)', linestyle='--')
plt.title('Autofunción correspondiente al QNM fundamental', fontsize=16)
plt.xlabel('r*', fontsize=15)
plt.ylabel('Ψ(r*)', fontsize=15)
plt.legend(fontsize=16) # Tamaño de fuente notablemente mayor para la leyenda
plt.grid(True)
plt.tick_params(axis='both', which='major', labelsize=13)
plt.tight_layout()
plt.show()