Author: Jonah Miller (jonah.maxwell.miller@gmail.com)
A Pseudospectral Elliptic Solver for Axisymmetric Problems Implemented in Python
Simply clone the repository and use
python setup.py install
Pyballd uses Chebyshev pseudospectral derivatives to attain very high accuracy with fairly low resolution. For example, if we numerically take second-order derivatives of this function:
and vary the number of points (or alternatively the maximum order of polynomial used for differentiation), we find that our error decays exponentially with the number of points. This is called "spectral" or "evanescent" convergence:
The appropriate domain for an axisymmetric problem is
where rh is some minimum radius. Infinite domains are difficult to handle. Therefore, we define
for most boundary situations or (following the work of [1])
when rh is a black hole event horizon. Then, following Boyd [2], we then define
for some characteristic length scale L so that X is defined on the domain [-1,1]. We perform our differentiation on X, which has no effect on the original PDE system except the introduction of Jacobian terms of the form
in a few places. Since one may want to assume additional (or different!) symmetry in the longitudinal direction, we do not impose any restriction there.
When x is defined as
the Jacobian for the coordinate transformation looks like
As we shall see, this spectral method can efficiently represent both exponential and algebraic decay with spectral accuracy.
The convergence of the scheme depends senitively on the characteristic length scale L. Generically it will be spectral but subgeometric, meaning it's not quite as fast as in the non-compact case.
Consider this function
which has this derivative
On the compact domain (on the equator), this function becomes
with derivative
We achieve best convergence for this function with L=1:
On the other hand, consider this function:
which has this derivative
On the compact domain (on the equator), this function becomes
with derivative
We achieve best convergence for this function with
as shown here:
In Pyballd, an elliptic system is defined via a residual. A residual
acts on a state vector u and its first and second derivatives in (in our case, axisymmetry) r and θ. If
then u(r,θ) is a solution to the PDE system.
An elliptic PDE is not well-posed without the addition of boundary conditions, which select for the particular solution. At infinity (X=1), we automatically demand that the solution must vanish. (In other words, we demand that all solutions are square-integrable.)
Dirichlet, Neumann, or Robin boundary conditions can be imposed on the inner radius (r_h), the position of minimum θ (often the axis of symmetry), and the position of maximum θ (often the axis of symmetry or the equator).
In the case of black-hole like compactifications, where
we automatically impose Von-Neumann boundary conditions at the inner boundary such that
which is a regularity condition we need to impose on solutions in these coordinates.
Pyballd simply requires that the user pass in functions that vanish when the residual and boundary conditions are satisfied. Along with information about the domain, such as the value of rh, this is sufficient information to construct a solution.
Poisson's equation is
or, in spherical coordinates it is
If we assume axisymmetry so that the azimuthal derivatives vanish and multiply both sides of the equation by the appropriate factors, we attain
For simplicity, we further restrict ourselves to the source-free case and attain
We would like to solve this problem using pyballd.
We begin by importing pyballd
and numpy
.
import pyballd
import numpy as np
pyballd
expects a residual, which defines the PDE system. The PDE is
satisfied when the residual vanishes. We define ours as
def residual(r,theta,u,d):
u = u[0]
out = (2*np.sin(theta)*r*d(u,1,0)
+ r*r*np.sin(theta)*d(u,2,0)
+ np.cos(theta)*d(u,0,1)
+ np.sin(theta)*d(u,0,2))
out = out.reshape(tuple([1]) + out.shape)
return out
Here d
is a derivative operator. So d(u,2,0)
corresponds to two
derivatives with respect to r of u, while d(u,0,2)
corresponds
to two derivatives with respect to θ.
The weird array reshaping is an artifact of the fact that pyballd
is
designed for vector and tensor equations, not just scalar
equations. Therefore, scalars must be treated as vectors of length 1.
pyballd
allso expects a boundary condition for the inner boundary. We
choose a Dirichlet boundary condition and define it as:
k = 4
a = 2
def bdry_X_inner(theta,u,d):
u = u[0]
out = u - a*np.cos(k*theta)
out = out.reshape(tuple([1]) + out.shape)
return out
This works a lot like residual
defined above. However, it will only
be evaluated at the inner boundary. The boundary condition is
satisfied when bdry_X_inner
vanishes.
Nonlinear elliptic systems are not necessarily unique. (And even linear ones may be difficult to uniquely solve numerically.) Therefore, we must feed the solver with an initial guess for the solution. We define ours as
def initial_guess(r,theta):
out = 1./r
out = out.reshape(tuple([1]) + out.shape)
return out
which is a simple square integrable function. It's not a solution that matches the boundary conditions. However, it should be sufficiently close to the true solution to allow for convergence.
Finally, we ask pyballd
to generate a solution by calling
pyballd.pde_solve_once
. Because it's a spectral method, you don't
need many nodes! For example:
SOLN,s = pyballd.pde_solve_once(residual,
r_h = 1.0,
order_X = 24,
order_theta = 6
theta_max = np.pi/2,
bdry_X_inner = bdry_X_inner,
initial_guess = initial_guess)
pyballd.pde_solve_once
returns the solution on the product grid of
colocation points and a discretization object s
, which can
differentiate a function defined on the colocation points, or
interpolate it to a uniform grid. For example:
SOLN = SOLN[0]
r = np.linspace(r_h,4,200)
theta = np.linspace(0,np.pi/2,200)
R,THETA = np.meshgrid(r,theta,indexing='ij')
interpolator = s.get_interpolator_of_r(SOLN)
soln_interp = interpolator(R,THETA)
x,z = R*np.sin(THETA), R*np.cos(THETA)
plt.pcolor(mx,mz,soln_interp)
plt.xlim(0,2)
plt.ylim(0,2)
will produce a beautiful plot of the solution interpolated to a finer grid, such as this one:
Both pde_solve_once
and discretization objects, called
PyballdStencil
s in the code, have a large number of options and
defaults. For a full description of of these, see the help strings
associated with them.
[1] Herderio, Radu, Runarrson. "Kerr black holes with Proca hair." Classical and Quantum Gravity 33-15 (2016).
[2] Boyd, John P. "Orthogonal rational functions on a semi-infinite interval." Journal of Computational Physics 70.1 (1987): 63-88.