Dear community, I got the bug!
The first definition I gave was indeed correct besides of a factor
cx=2pi/N, as Steven pointed out. Forget about the second definition.
The bug was in the sign of the solution of the dif eq that I wrote in the
PD. This sign appears because of the nabla operator: -ik, and has to be in
agreement with the fftw definition of the FFT and not with your own
one!!!!
Thank you everybody!
Cheers!
Correct definition:
>PD:
>The solution of the dif eq in 2d is then:
>
> for (int i=0;i<Nx;i++){
> for (int j=0;j<Ny;j++){
> re(vx[j+Ny*i]) = c* kx[i]/(kx[i]*kx[i]+ky[j]*ky[j])
Sorry, I misread your code, I thought your first loop went from i <
Nx/2 when in fact you used i < Nx.
Reply by paquillo●March 20, 20062006-03-20
Thank you for your comment Steven.
however let me say you that I am still confused.
1)"kx[Nx/2] and ky[Ny/2] are not defined (for even Nx and Ny)", why, what
can I do?
2)Now, d(x) is a Nx*Ny=Nx*Nx array: d(0) ... d(Nx*Nx-1)
The corresponding kx component should then go from kx(0) till kx(Nx-1), so
I do not understand why you tell me to put <= in the "first loop".
3)I imposed Nx=Ny (always even) with an equidistant grid, and thus kx=ky.
The normalization 2*pi/(Nx*dx) should not affect the qualitative result.
Besides that I thought this normalization is unnecessary since it is
carried in the exponent in the definition of the fftw.
4)I fully agree with you in the final comment. I took this into account.
For the first definition of the k_vector I imposed v=0 if k=0. For the
second definition this does not happen (k=!0).
Cordially,
paco
Reply by ●March 19, 20062006-03-19
paquillo wrote:
> Now I would expect the following defintion of the k_vector to be correct:
> for (int i=0 ; i<Nx;i++){
> kx[i] = i;
> }
> for (int i=1 ; i<=Nx/2;i++){
> kx[Nx-i] = -kx[i];
> }
This looks mostly fine at first glance, except kx[Nx/2] and ky[Ny/2]
are not defined (for even Nx and Ny): you need a "<=" on your first
loop. Note also that you should really use 2*pi*i/(Nx*dx) where dx is
your discretization in the x direction, and similarly for y; this
matters if your computational cell has a different size in the x and y
directions since otherwise your kx and ky won't have compatible units.
Also, maybe you forgot to special-case the (0,0) point when you are
solving the PDE, since if you are not careful you'll end up dividing by
zero there.
Cordially,
Steven G. Johnson
Reply by paquillo●March 19, 20062006-03-19
Hello everybody!
I am having problems with the defintion of the k vector trying to solve a
2 and 3 dimensional differential equation in Fourier space, namely:
div v_vector(x)=c*d(x).
This can be written in Fourier space like:
v_vector(k)=i**c*d(k)*k_vector/k^2.
I am doing the complex to complex fftw transform of d(x) adding zero
imaginary part to d(x).
Now I would expect the following defintion of the k_vector to be correct:
for (int i=0 ; i<Nx;i++){
kx[i] = i;
}
for (int i=0 ; i<Ny;i++){
ky[i] = i;
}
for (int i=1 ; i<=Nx/2;i++){
kx[Nx-i] = -kx[i];
}
for (int i=1 ; i<=Ny/2;i++){
ky[Ny-i] = -ky[i];
}
This means, that I have first the 0 frequency, then from 1 till the
nyquist frequency, N/2, and then from -N/2 til -1.
However, this defintions give wrong results (unreasonable results),
instead the following definition seems to work (very reasonable results)
and I do not understand why , since it should be completely wrong:
for (int i=0 ; i<Nx;i++){
kx[i] = i+1;
}
for (int i=0 ; i<Ny;i++){
ky[i] = i+1;
}
for (int i=0 ; i<Nx/2;i++){
kx[Nx-1-i] = kx[i];
}
for (int i=0 ; i<Ny/2;i++){
ky[Ny-1-i] = ky[i];
}
Does somebody know which is the correct definition of the wave vector k in
Fourier space?
Thank you in advance!
PD:
The solution of the dif eq in 2d is then:
for (int i=0;i<Nx;i++){
for (int j=0;j<Ny;j++){
re(vx[j+Ny*i]) = c* kx[i]/(kx[i]*kx[i]+ky[j]*ky[j]) im(d[j+Ny*i]);
im(vx[j+Ny*i]) = c* kx[i]/(kx[i]*kx[i]+ky[j]*ky[j])-re(d[j+Ny*i]);
re(vy[j+Ny*i]) = c* ky[j]/(kx[i]*kx[i]+ky[j]*ky[j]) im(d[j+Ny*i]);
im(vy[j+Ny*i]) = c* ky[j]/(kx[i]*kx[i]+ky[j]*ky[j])-re(d[j+Ny*i]);
}
}