Using spectral methods to solve the Helmholtz problem
This code is a tutorial code for solving the Helmholtz (Poisson) problem in a doubly periodic domain. Since the diffusivity is constant, once we compute the Fourier transforms, solving for is trivial.
Define the various grid parameters.
Nx = 256; % uonal resolution Ny = 256; % meridional resolution Lx=3e3; Ly=3e3; dx=Lx/Nx; dy=Ly/Ny; x=dx*(1:Nx)'; y=dy*(1:Ny)'; dx2o2 = dx^2/2; dy2o2 = dy^2/2; [xx,yy] = meshgrid(x,y);
FFT related material
For this problem we will be solving the PDE using spectral methods. We transform the problem using Fourier transforms so that the derivatives become multiplication.
Build the wavenumbers and create the variable diffusivity.
dk = 2*pi/Lx; % Nyquist wavenumbers dl = 2*pi/Ly; k=[0:Nx/2-1 Nx/2 -Nx/2+1:-1]'*dk; l=[0:Ny/2-1 Ny/2 -Ny/2+1:-1]'*dl; [kk,ll]=meshgrid(k,l); l2 = l.^2; % 1D meridional wavenumber arrays for splitting up the 1D problems il = 1i*l; k2 = k.^2; % 1D uonal wavenumber arrays for differentiating within a sub problem ik = 1i*k; ll2 = ll.^2; % 2D meridional wavenumber arrays for differentiation purposes ill = 1i*ll; kk2 = kk.^2; % 2D uonal wavenumber arrays for general differentiation purposes ikk = 1i*kk;
The diffusivity is constant here. In another script we will explore what happens when it is non-constant.
kappa = 100;
Solve the problem via multiplication in spectral space
We solve for u by taking the Fourier transform of our RHS and multiplying it by . We then invert the transform to return to physical space.
RHS = sin(6*pi*xx/Lx).*sin(4*pi*yy/Ly); ftRHS = fft2(RHS); multfac = 1./(kappa*(kk2+ll2)); ftu = multfac.*ftRHS; ftu(1,1)=ftRHS(1,1); u = real(ifft2(ftu));
Plot the Results