I want to solve Poisson equation in the ellipse via (pseudo)-spectral methods.

I was thinking of using Dedalus via Pycall or something like that. But are there any good alternatives in Julia? Most things that I’ve found are for periodic boundary conditions or 1D but I really need the curvilinear option.

The ApproxFun.jl package is based on banded-matrix spectral techniques similar to those in Dedalus, though it doesn’t support quite the same range of domains AFAIK.

Here is an example of a Helmholtz equation on a disk with the MultivariateOrthogonalPolynomials package (by the same authors as ApproxFun), that you could probably adapt to Poisson on an elliptical disk: MultivariateOrthogonalPolynomials.jl/diskhelmholtz.jl at master · JuliaApproximation/MultivariateOrthogonalPolynomials.jl · GitHub

That being said, I don’t think there is currently anything (in any language) that matches Dedalus for scaling spectral methods to huge parallel calculations.

We’re planning to work on the disk over the next couple months. There will be some places where we will have more features than Dedalus, eg variable coefficient Helmholtz, but if you want to do time-evolution with parallelisation we won’t be anywhere close to Dedalus (as it’s not my main research interest).

Note for time-evolution using splitting methods the biggest cost is the transforms. FastTransforms.jl has a zernike polynomial transform so using that it wouldn’t be too hard to implement a simple splitting method.

I just manually discretize the spatial region with Chebyshev or Fourier collocation points.

Using collocation on a disk with radial coordinates is problematic as it introduces an artificial singularity at the origin and will be badly conditioned. But if you your right hand side is nice it will work.

Better would be to follow Wilber-Townsend-Wright which incorporates symmetries and uses coefficient space methods. This still has weak singularities at the origin but for Poisson that is fine. Unfortunately I don’t know of a Julia implementation though the transform part is in FastTransforms.jl.

Zernike polynomials fully diagonalise the laplacian with Dirichlet conditions and so the solve is very well-conditioned with no artificial singularities.