# build a simple incompressible solver with spectral method

I am trying to build a simple solver for the incompressible fluid in a periodic box to learn the spectral method. I am following the textbook (Peyret R. Spectral methods for incompressible viscous flow[M]. Springer Science & Business Media, 2013.) the author presents a Fourier method for the calculation of spatially periodic solutions. The first step is applying Fourier transform to momentum and continuity equation:

($A$ is nonlinear term A(U,U), $f$ is source term)

Then apply divergence operator in the Fourier space to momentum equation, using continuity equation in Fourier space to eliminate $p$ in momentum equation:

It looks like it eliminates the necessity of projecting and solving pressure for each time step because the transformed equation satisfied the continuity automatically. However, when I tried to implement the code in Matlab for the 2d periodic domain, I didn't get a reasonable result and I have no idea about how to fix it.

My question here is, did I miss something nontrivial for implementing this method? Or just some bugs in my code?

(this is only part of the code to keep brevity)

%% Construct k
k1 = 1i*[0:N/2-1 0 -N/2+1:-1];
k2 = -([0:N/2 -N/2+1:-1]).^2;

% expand to 2d
k1 = repmat(k1,N,1);
k2 = repmat(k2,N,1);

%% Solve PDE
for n = 1:plotgap

t = t+dt;

%% ifft velocity field
u_r = real(ifft2(u));v_r = real(ifft2(v));

ux_r = real(ifft2(k1.*u)); uy_r = real(ifft2(k1'.*u));
vx_r = real(ifft2(k1.*v)); vy_r = real(ifft2(k1'.*v));

A1 = fft2(u_r.*ux_r + v_r.*uy_r);
A2 = fft2(u_r.*vx_r + v_r.*vy_r);

% Mass conservation
C1 = (k1.*A1 + k1'.*A2)./(k2 + k2') .*k1;
C2 = (k1.*A1 + k1'.*A2)./(k2 + k2') .*k1';
C1(1,1) = 0;C2(1,1) = 0;

% Diffusion
D1 = nu*(k2 + k2') .* u;
D2 = nu*(k2 + k2') .* v;

%% Construct operator
operator1 = -A1 + D1 + C1;
operator2 = -A2 + D2 + C2;

%% Update velocity

if ~exist('tmp1','var'), tmp1 = operator1; end
if ~exist('tmpp1','var'), tmpp1 = operator1; end
if ~exist('tmppp1','var'), tmppp1 = operator1; end
if ~exist('tmp2','var'), tmp2 = operator2; end
if ~exist('tmpp2','var'), tmpp2 = operator2; end
if ~exist('tmppp2','var'), tmppp2 = operator2; end

u = u + dt*(55/24*operator1 - 59/24*tmp1 + 37/24*tmpp1 - 3/8*tmppp1);  % 4th -order Adam - bashforth
tmppp1 = tmpp1;tmpp1 = tmp1;tmp1 = operator1; % update f values
v = v + dt*(55/24*operator2 - 59/24*tmp2 + 37/24*tmpp2 - 3/8*tmppp2);
tmppp2 = tmpp2;tmpp2 = tmp2;tmp2 = operator2; % update f values

end


1. The Nyquist frequency N/2 should be present in both k2 and k1 (or absent from both)
2. k1 = 1i*[0:N/2 -N/2+1:-1]; [KX,KY] = meshgrid (k1); K2=KX.^2+KY.^2 would be a much cleaner way to define all your wavenumber arrays: the wavenumber along x, KX, along y KY, and its square K2. This is not only a question of taste, if would clean your code a lot by eliminating the need to use transposes.
3. Speaking of transposes, ' is a complex conjugate transpose whereas .' is a regular transpose. Computing k1' does not do what you think it does in your code.