# Integrating the 2d vorticity equation on periodic boundaries

This question is a follow-up of https://stackoverflow.com/questions/44718160/solve-a-linear-system-for-fft-coefficients

At some time (kt), the FT of the vorticity (omega) satisfies:

(1-nu*L*dt/2)fftomega(kt+1)=(1+nu*L*dt/2)fftomega(kt)-fft_J*dt


where fftomega is the 2d fft of omega, J(psi,omega) is the Jacobian (see below) and the streamfunction del^2(psi)=omega. L is an array of coefficients that when multiplied by say fftomega return the fft of the laplacian of omega. For a 6x6 L is

0 -1 -4 -9 -4 -1

-1 -2 -5 -10 -5 -2

-4 -5 -8 -13 -8 -5

-9 -10 -13 -18 -13 -10

-4 -5 -8 -13 -8 -5

-1 -2 -5 -10 -5 -2

The fft of the Jacobian is computed as follows (using matlab/octave notation for brevity):

Dx_psi=real(ifft2(Dx.*fftpsi));                    %v
Dy_psi=real(ifft2(Dy.*fftpsi));                    %-u
Dx_omega=real(ifft2(Dx.*fftomega));
Dy_omega=real(ifft2(Dy.*fftomega));
Jacobian = - Dy_psi.*Dx_omega + Dx_psi.*Dy_omega;

fft_J = fft2(Jacobian);


In practice, the Jacobian is integrated with a predictor corrector scheme due to Gadzag (the -fft_j*dt in the equation above is schematic). If arrays such as fftomega have the same dimension as omega, then the equation above can be solved as a system of linear equations, with LAPACK for instance. However since the physical variables are real, procedures such as FFTW allow for output fft's to be non-square (r2c,c2r). In this case it seems like the only way to integrate the equation is to expand these fft's back out to square, giving up the advantage afforded by the r2c,c2r scheme. Alternatively, one could make the physical variables complex with null imaginary part, and use the c2c scheme.

Has anyone found an alternative that makes it possible to use r2c,c2c?

• Note: the question is a follow-up of stackoverflow.com/questions/44718160/… Jun 23, 2017 at 15:33
• I still don't understand the difference from the last case I answered. There is no matrix here. I don't see what you wan't to do with LAPACK. fftJ(i) is again just a number for each wavenumber and all wavenumbers are independent. Please clarify. I asked for the details and we really need them. Jun 23, 2017 at 15:35
• You were talking about some structure of some matrix last time. Which matrix? How does the structure look like? Maybe it is the fftJ and I understand it wrong, but then a matrix should be somehow applied on the vector of unknowns. How? Jun 23, 2017 at 15:37
• @VladimirF I added a description of what fft_J is. This is an aliased, version, in practice I use an anti-aliasing algorithm Jun 23, 2017 at 16:12
• It still looks like a simple number for each Fourier mode, no matrix, no linear system, no LAPACK. Jun 23, 2017 at 16:57