build a simple incompressible solver with spectral method

Question

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

    %% Time advance
    t = t+dt;

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

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

    % Advection
    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

Answer 1

Your algebra is fine. Periodic boxes are particularly nice for imposing incompressibility: as you found out, you can easily get rid of pressure. [This intimately has to do with the fact that differential operators are diagonal on the Fourier basis, so that solving an otherwise horrendous Laplace equation for pressure is straightforward here.]

A couple remarks on the implementation:

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.
4. Your time-stepping scheme seems to be a 4th order Runge-Kutta. I would start with a naive Euler, or Backward Euler, debug that, and then ramp up the complexity of the time-stepping scheme.

I find this document to be a goldmine as an initiation to spectral methods and computational fluid dynamics. All the fundamentals are there. If you take the time to read it, understand it, and write your own programs following the guidance it provides, you will have a very solid foundation for writing spectral codes.