Speed of ApproxFun.jl solution to Screened Poisson in cylindrical coordinates

Specific Domains Modelling & Simulations
1

The examples of ApproxFun.jl examples solutions are fascinating and powerful, and I’d love to figure out the best performance for solving the homogeneous screened Poisson equation but in cylindrical coordinates (with circular symmetry) and with a mix of Dirichlet and homogeneous Neumann boundaries.

\Delta u(r,z) - \lambda^2 u(r,z) = 0 \\ u, \lambda \in \mathbb{R}, \quad r \in [a,1], \quad 0<a<1, \quad z \in [0, 1] \\ \frac{\partial}{\partial z}u(r,z)|_{z=0} = 0 \quad u(r,1) = k_1 \\ u(a,z) = k_2 \quad u(1,z) = k_3.

Here with the particular cylindrical Laplacian and circular symmetry (\tfrac{\partial}{\partial \theta} = 0):

\Delta u(r,z) = \frac{1}{r} \frac{\partial}{\partial r}\left( r \frac{\partial}{\partial r}u(r,z) \right) + \frac{\partial^2}{\partial z^2}u(r,z) \\ = \frac{1}{r} \frac{\partial}{\partial r}u(r,z) + \frac{\partial^2}{\partial r^2}u(r,z) + \frac{\partial^2}{\partial z^2}u(r,z) .

The ApproxFun approach below seems to solve this problem, albeit a bit slowly. It is slower when setting constant Dirichlet boundaries than when using the \exp(r)\cos(z) expressions. I’ve reduced the data using resizedata!and reduced the tolerance with tolerance=1E-5.

λ = 10
a = 0.5
k_1 = 1.2
k_2 = 0.2
k_3 = 1.0
dr = a .. 1.0
dz = 0.0 .. 1.0
d = dr × dz
r, z = Fun(∂(d))
r, z = components(r), components(z)
g = [-exp(r[1])*sin(0); exp(1)*cos(z[2]);
      exp(r[3])*cos(1);  exp(a)*cos(z[4])]
# g2 = [0; k_3; k_1; k_2]
Δ = Laplacian(d)
Dr = Derivative(d, [1, 0])
B = [I⊗lneumann(dz); rdirichlet(dr)⊗I; I⊗rdirichlet(dz); ldirichlet(dr)⊗I]
invfun = Fun((r, z) -> 1 / r, d)
L = [B; invfun * Dr + Δ - λ^2 * I]
Q = qr(L)
ApproxFunBase.resizedata!(Q,:,10);
u = \(Q, [g; 0]; tolerance=1E-5)

What else can I do to speed up the solution?

PS: Thank you for making these tools in the first place.

2

You can put this in a function but I am not sure you’ll gain much. Is depends on how much time is spent in u = \(Q, [g; 0]; tolerance=1E-5).

I dont use AF for its speed but for its precision. It is fantastic to debug less precise code, or to compute the spectrum of operators with very high precision, etc.

3

Thank you - putting the step u = \(Q, ...) in a function does speed up many of the solutions.

Still, some choices of constant Dirichlet boundary conditions don’t finish in the several hours my patience lasts, so I hope there’s still some way to do with ApproxFun.