Losing radial symmetry in solving Poisson equation with FFTW in C

Question

I have a problem about using FFTW in order to solve Poisson equation. The equation I am trying to solve is the following:

The problem is that I get a solution which is no longer radially symmetric, in particular symmetry is lost along the axes as you can see from the output plot:

Here's the code:

    double dx = 25. / ( N*1. ), dy = 25. / ( N*1. ), dkx = 2*pi/(N*dx), dky = 2*pi/(N*dy);//define steps in real and Fourier space
    double Lx = N*dx, Ly = N*dy, Lkx = N*dkx, Lky = N*dky;//define box sides in real and Fourier space
    int i, j, Nh = N/2 + 1;
    double *inr, *in2r, *in2d;
    fftw_complex *outr;
    fftw_complex *outd;
    fftw_plan rplan_forward;
    fftw_plan dplan_backward;
    fftw_plan rplan_backward;

    for ( i = 0; i < N; i++ ) 
    { 
        x[i] = dx*i - Lx/2.;
        kx[i] = dkx*i - Lkx/2.;
        y[i] = dy*i - Ly/2.;
        ky[i] = dky*i - Lky/2.;
    }

    outd = (fftw_complex*)fftw_malloc ( sizeof ( fftw_complex ) * N * Nh );

    inr = ( double * ) malloc ( sizeof ( double ) * N * N );


    for ( i = 0; i < N; i++ ) //Initialize ine to rho
    {
        for ( j = 0; j < N; j++ )
        {
            inr[i*N+j] = rho[i*N+j] * pow( -1, i + j );                            
        }
    }

    outr = (fftw_complex*)fftw_malloc ( sizeof ( fftw_complex ) * N * Nh );

    rplan_forward = fftw_plan_dft_r2c_2d ( N, N, inr, outr, FFTW_ESTIMATE );

    fftw_execute ( rplan_forward ); //fft rho

    for ( i = 0; i < N; i++ ) //Evaluate d in Fourier Space according to Poisson eq
    {
        for ( j = 0; j < Nh; j++ ) 
        {
            outd[i*Nh+j][0] = g * outr[i*Nh+j][0] / ( kx[i]*kx[i] + ky[j]*ky[j] ); 
            outd[i*Nh+j][1] = g * outr[i*Nh+j][1] / ( kx[i]*kx[i] + ky[j]*ky[j] ); 
        }                      
    }

    outd[Nh*N/2+N/2][0] = 0.; //Avoid singularity at k = 0
    outd[Nh*N/2+N/2][1] = 0.;

    in2d = ( double * ) malloc ( sizeof ( double ) * N * N );

    dplan_backward = fftw_plan_dft_c2r_2d ( N, N, outd, in2d, FFTW_ESTIMATE );

    fftw_execute ( dplan_backward ); //Antitrasform d

    in2r = ( double * ) malloc ( sizeof ( double ) * N * N );

    rplan_backward = fftw_plan_dft_c2r_2d ( N, N, outr, in2r, FFTW_ESTIMATE );

    fftw_execute ( rplan_backward ); //Antitrasform rho

    for( i = 0; i < N; i++ ) 
    {
        for ( j = 0; j < N; j++ ) 
        {
            printf( " %le %le %le \n", x[i], y[j], in2d[i*N+j] * pow( -1, i + j ) / ( N*N*1. ) - in2d[0] / ( N*N*1. ) ); //print d in such a way that it is zero at the boundaries
        }

    }

    fftw_destroy_plan ( rplan_forward );
    fftw_destroy_plan ( dplan_backward );
    fftw_destroy_plan ( rplan_backward );
    fftw_free ( outd );
    fftw_free ( inr );
    fftw_free ( outr ); 
    fftw_free ( in2d ); 
    fftw_free ( in2r); 
    fftw_cleanup();

Notice that I would like to have a solution which goes to zero at the boundaries.

Could you please tell me where is the problem? I have the suspect that it is just a numerical issue but I am not sure about it.

Here is a zoom on the boundaries ( in particular on the region -12.5 < y < -11.5 ).

Answer 1

The loss of radial symmetry is related to the shape of the domain : it is a square, not a circle. Furthermore, since the dicrete Fourier transform is applied, the problem that is studied is the Poisson's equation on a square domain using periodic boundary conditions. Hence, the temperature near the boundaries may not decrease to zero.

This is the reason why the zero frequency, that is the average of the source term, is not used. If it were not the case, the stationnary heat equation could not be solved. Indeed, due to the periodic boundary conditions, the heat flux getting out on one side is getting in on the opposite side. As when the domain is perfectly insulated, the temperature would go up and never reach a stationnary solution.

To apply zero Dirichlet boundary solution, one of the sine transforms must be used instead of the Fourier transform. See for instance:

http://www.ds.ewi.tudelft.nl/vakken/in4049TU/sheets/lecture10.pdf

Here is an answer (of mine!) featuring a piece of code solving 2D Poisson's equation, zero Dirichlet boundary condition, by mean of sine transform of FFTW. It's c++, but translating it to c is straightforward.

fftw3 for poisson with dirichlet boundary condition for all side of computational domain

To recover something close to a radial symmetry, the support of the source must be far from the edges of the domain. The best solution might be to switch to a 1D radial equation, if you know that the domain, boundary conditions and the source term respect the radial symmetry.

Answer 2

I think one issue is here:

outd[Nh*N/2+N/2][0] = 0.; //Avoid singularity at k = 0

k=0 is at (0,0), not (N/2,N/2). This is how the DFT is defined. k runs from 0 to N-1. You are introducing a phase shift in the spatial domain. I don’t know how this affects your later computations. You are trying to center the spectrum by multiplying the input with pow(-1,i). This trick works only for even N!

If is better to compute your kx and ky properly:

for ( i = 0; i < N; i++ ) {
   kx[i] = ( i > N/2 ? i-N : i ) * 2*pi/N;
}

This will work for any N. Also note that my definition of k above corresponds to the radial frequency (ω), not the integer k of the DFT. It also differs from yours, I think (if so, yours is incorrect!).

The other issue is that the DFT is not the FT. The DFT is computed on a bounded spatial domain, which is presumed periodic. This domain is a square. It is hard to do rotationally symmetric things on a square grid. The Nyquist frequency (max possible frequency) is sqrt(2) times as high along the diagonal compared to the axes. And the period of the pattern is also sqrt(2) times as high along the diagonal. Thus, you get more aliasing along the axes, and your pattern is more likely to overlap itself along the axes.

If the above is causing the effect you’re seeing, try increasing the size of your spatial domain, and/or the sampling density.