3D Fourier-spectral solver for Poisson's equation with mixed boundary conditions is only first-order convergent.

[ali-ramadhan]

The 3D solver with mixed boundary conditions (periodic in the horizontal with DFTs, Neumann in the vertical with DCTs) works now as tested against an analytic solution but for some reason once I switched to DCTs it's only first-order convergent. The relevant function is solve_poisson_3d_mbc.

[glwagner]

Could this mean that the boundary condition is incorrect? I believe we expect first-order convergence in that case --- for example, projecting a step function onto a Fourier series.

[ali-ramadhan]

Hmmm that might be it.

I tried to construct an analytic solution that obeys all the boundary conditions: using the source term f(x,y,z) = 4(x²+y²-1-π²/H²) cos(2πz/Lz) exp(-(x²+y²)) the analytic solution is ϕ(x,y,z) = cos(2πz/H) exp(-(x²+y²)) which does satisfy ∂ϕ/∂z = 0 at z=0 and z=H, and is "numerically periodic" as the exponential should decay to < machine epsilon at the boundaries in a large enough domain, but maybe I did something wrong in the process.

I'll try testing with a different analytic solution. I wanted to avoid a function that is just sines and cosines as that would give a spectral solution that is accurate to within machine precision, which doesn't help with testing convergence properties.

[glwagner]

Is it true that the different cosine transforms differ in their definition of the collocation points?

Where is the pressure stored in z? At cell centers? So your grid extends from z=-dz/2 to z=-H+dz/2? Just from glancing at the FFTW documentation those would seem the right grid points for DCT-II.

You can have more than one test! I think it makes sense to have an easy test with sinusoids to catch obvious bugs. Then if you are ambitious you can also implement a Gaussian test.

[ali-ramadhan]

Yes the pressure is stored as the cell centers so the DCT-II seems to be the right way to go, which is nice as FFTW.dct() uses the DCT-II. Then I assumed DCT-III is the correct inverse transform, which is what FFTW.idct() uses, but DCT-III is even around j=0 and odd around j=n while DCT-IV is even around j=-0.5 and odd around j=n-0.5, so maybe we should be using DCT-II for the forward and DCT-IV for the

I'm also trying out a spectral solve in the horizontal followed by an explicit tridiagonal solve in the vertical to see if I can get it to work.

But yeah maybe that's why the Gaussian test didn't work as it's even around j=0.

The Gaussian test worked in 2D so it would make sense that the DCT collocation points are now resulting in first-order convergence.

[glwagner]

I wrote a tridiagonal solver here:

https://github.com/LilithCFD/Lilith.jl/blob/master/src/solvers.jl

and a test:

https://github.com/LilithCFD/Lilith.jl/blob/master/test/runtests.jl

[glwagner]

I'm not so sure about DCT-IV being correct --- the j's for the inverse transform refer to the spectral-space wavenumbers, not the physical collocation points.

The part of the solution corresponding to j=0 (in FFTW's notation) should be the barotropic mode; aka the vertical average. That doesn't seem to be the case for DCT-IV?

[ali-ramadhan]

Ah you're right, the DCT-III's j=0 refers to the wavenumbers so DCT-III is correct. I suppose this makes sense as the DCT-II and III always seem to go together.

Thanks for the Julia tridiagonal solver! Hope you don't mind if I use it here.

Not sure if it applies to this problem but the tridiagonal matrix algorithm

is not stable in general, but is so in several special cases, such as when the matrix is diagonally dominant (either by rows or columns) or symmetric positive definite.

[glwagner]

Yes, it is definitely ok for you to re-use, copy-paste, or declare-completely-wrong any code I have!

Another tridiagonal implementation is used by Diablo. Maybe that will be useful.

[ali-ramadhan]

Closing as this particular solver is not used anymore (it was wrong).