[jollage]

Hi,

I would like to do dealiasing for 2D data.

I have two 2D arrays of size, say N=5. I want to multiply them, but before that I would like to do the so-called 3/2 zero padding for dealiasing (following Orszag 1971). This is the famous pseudospectral method for treating nonlinear terms in direct numerical simulations. The code with comments below shows how I did it, I would like to know if I'm doing it correctly. Thanks.

I work in Matlab. Here is my code

N=5;
a_p=rand(N,N);
b_p=rand(N,N);

% go to spectral space
a_s=fft2(a_p)/N/N;
b_s=fft2(b_p)/N/N;

% empty columns and rows are inserted for padding zero in the spectrum
aa_aux=[a_s(:,1: (N-1)/2+1) zeros(N,(N-1)/2) a_s(:,(N-1)/2+1+1:N)];
a_s=[aa_aux(1: (N-1)/2+1,: );zeros((N-1)/2,N+(N-1)/2);aa_aux((N-1)/2+1+1:N,: )];

% empty columns and rows are inserted for padding zero in the spectrum
bb_aux=[b_s(:,1: (N-1)/2+1) zeros(N,(N-1)/2) b_s(:,(N-1)/2+1+1:N)];
b_s=[bb_aux(1: (N-1)/2+1,: );zeros((N-1)/2,N+(N-1)/2);bb_aux((N-1)/2+1+1:N,: )];

% go back to physical space
a_p2=ifft2(a_s)*(3*N-1)/2*(3*N-1)/2;
b_p2=ifft2(b_s)*(3*N-1)/2*(3*N-1)/2;

c_p=a_p2.*b_p2;

c_p=c_p(1:N,1:N);

But when I compare the "dealiased" c_p with the straightforward result a_p.*b_p, I can't link them in any sense except that the first mode is the same. Any comments on this? Or did I do something wrong? Thanks a lot in advance.

[kaya]

Two suggestions :

You should try to do this on 1D first to get the idea, if you do this in 1D I am sure that you'll spot your mistakes. If you want to learn you should do this. If you don't care and just want to make it work, I have some matlab routines that I can pass.

The order should be like

%take ffts
ah=fft(a)
bh=fft(b)
%zero pad
ah0=zeropad(ah)
bh0=zeropad(bh)
%get back to real space
a0=ifft(ah0)
b0=ifft(bh0)
%do the operation
c0=a0*b0
%go back to fourier space
ch0=fft(c0)
%take the contaminated waves out
ch=trim(ch0)
%normalize
ch=normalize(ch)
%get back to real space
c=ifft(ch) %dealiased product


Don't give random values, use sin cos with known amplitudes that you specify i.e. 2*sin(2x)+3*sin(4x)+...sin(mx) and make sure that m is equal to N/2 which will be your nyquist freq. == in other words, the highest frequency that you can possibly resolve on that grid.



Here some more lines since I am not sure if you know what you're doing.
Your padding in fourier space is nothing more than interpolation. You get back to the real space with a finer resolved wave on a grid of 3*N/2 and do your multiplication. You like to take your multiplication on an interpolated (finer) grid because your previous grid with N points was barely enough to resolve your high frequency waves, and the multiplication operation that you're doing will generate waves at even higher frequencies. From fourier space perspective: it is like if you don't have a place to put these waves that are about to be produced, they will contaminate the waves that you could resolve and the solution will become aliased.


also you need to know what matlab is doing when you call fft and ifft.
This one is a bit tricky. I haven't checkd your code one by one and it is usually try-and-see-if-it-works like thingy since its hard to spot the bugs by eye (at least for me). Here is the thing: when you take an fft(ones(4,1)) you will see your first amplitude will be 4. If you do the same operation with fft(ones(6,1)), then your first amplitude will be 6. So in fourier space the same field will be represented with different values just because one of them is on a finer grid. Matlab knows this and normalizes when you take ifft but since you change the grid in fourier space, it won't be able to correctly normalize. Which is the reason why you have to do this normalization book keeping well since you will be padding-unpadding (changing the grid// interpolating) in fourier space.

[jollage]

Originally Posted by kaya

Two suggestions :

You should try to do this on 1D first to get the idea, if you do this in 1D I am sure that you'll spot your mistakes. If you want to learn you should do this. If you don't care and just want to make it work, I have some matlab routines that I can pass.

The order should be like

%take ffts
ah=fft(a)
bh=fft(b)
%zero pad
ah0=zeropad(ah)
bh0=zeropad(bh)
%get back to real space
a0=ifft(ah0)
b0=ifft(bh0)
%do the operation
c0=a0*b0
%go back to fourier space
ch0=fft(c0)
%take the contaminated waves out
ch=trim(ch0)
%normalize
ch=normalize(ch)
%get back to real space
c=ifft(ch) %dealiased product


Don't give random values, use sin cos with known amplitudes that you specify i.e. 2*sin(2x)+3*sin(4x)+...sin(mx) and make sure that m is equal to N/2 which will be your nyquist freq. == in other words, the highest frequency that you can possibly resolve on that grid.



Here some more lines since I am not sure if you know what you're doing.
Your padding in fourier space is nothing more than interpolation. You get back to the real space with a finer resolved wave on a grid of 3*N/2 and do your multiplication. You like to take your multiplication on an interpolated (finer) grid because your previous grid with N points was barely enough to resolve your high frequency waves, and the multiplication operation that you're doing will generate waves at even higher frequencies. From fourier space perspective: it is like if you don't have a place to put these waves that are about to be produced, they will contaminate the waves that you could resolve and the solution will become aliased.


also you need to know what matlab is doing when you call fft and ifft.
This one is a bit tricky. I haven't checkd your code one by one and it is usually try-and-see-if-it-works like thingy since its hard to spot the bugs by eye (at least for me). Here is the thing: when you take an fft(ones(4,1)) you will see your first amplitude will be 4. If you do the same operation with fft(ones(6,1)), then your first amplitude will be 6. So in fourier space the same field will be represented with different values just because one of them is on a finer grid. Matlab knows this and normalizes when you take ifft but since you change the grid in fourier space, it won't be able to correctly normalize. Which is the reason why you have to do this normalization book keeping well since you will be padding-unpadding (changing the grid// interpolating) in fourier space.

Hi Kaya,

Thanks for your reply. I did 1d dealiasing first; for 1d, it's straightforward. I want to implement the 2D dealiasing. Your example is 1D dealiasing, could you pass me the 2D dealiasing routine you mentioned? Thanks.

I agree that I should have the additional steps of trimming the elongated spectrum back to the original one in the Fourier space.

I understand your concern about normalization. What I did just after the routines fft2 and ifft2 is for the normalization. Thanks.

[kaya]

c_p=c_p(1:N,1:N); %you want to do this in fourier space not in real

function [ ch ] = aap( ah,bh,deall )



if deall == 1
[ny, nz, nx] = size(ah); 

m  = nz*3/2;
mm = nx*3/2;

ahp = [ah(:,1:nz/2,:) zeros(ny,(m-nz),nx) ah(:,nz/2+1:nz,:)];
bhp = [bh(:,1:nz/2,:) zeros(ny,(m-nz),nx) bh(:,nz/2+1:nz,:)];

ahhp = shiftdim(ahp,1);
bhhp = shiftdim(bhp,1);

ahpn = [ahhp(:,1:nx/2,:) zeros(m,(mm-nx),ny) ahhp(:,nx/2+1:nx,:)]; % pad uhat with zeros 
bhpn = [bhhp(:,1:nx/2,:) zeros(m,(mm-nx),ny) bhhp(:,nx/2+1:nx,:)]; % pad uhat with zeros 
ahpn = shiftdim(ahpn,2);
bhpn = shiftdim(bhpn,2);

ap   =ifft(ifft(ahpn,[],3),[],2); 
bp   =ifft(ifft(bhpn,[],3),[],2); 

conv  = ap.*bp;
convh = fft(fft(conv,[],3),[],2);

ch = 1.5*[convh(:,1:nz/2,:) convh(:,m-nz/2+1:m,:)]; % extract F-coefficients 
ch = shiftdim(ch,1);
ch = 1.5*[ch(:,1:nx/2,:) ch(:,mm-nx/2+1:mm,:)]; % extract F-coefficients 
ch = shiftdim(ch,2);

elseif deall == 0 
a = real(ifft(ifft(ah,[],2),[],3));
b = real(ifft(ifft(bh,[],2),[],3));
c = a.*b;
ch = fft(fft(c,[],2),[],3);
end
end