Потеря радиальной симметрии при решении уравнения Пуассона с FFTW в C

Question

У меня проблема с использованием FFTW для решения уравнения Пуассона.Уравнение, которое я пытаюсь решить, выглядит следующим образом:

Проблема в том, что я получаю решение, которое больше не является радиально-симметричным, в частности симметриятеряется по осям, как вы можете видеть из выходного графика:

Вот код:

    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();

Обратите внимание, что я хотел бы получитьрешение, которое идет к нулю на границах.

Подскажите, пожалуйста, в чем проблема?У меня есть подозрение, что это просто численная проблема, но я не уверен в этом.

Здесь приведено увеличение границ (в частности, в области -12,5

Answer 1

Потеря радиальной симметрии связана с формой области: это квадрат, а не круг .Кроме того, поскольку применяется дискретное преобразование Фурье, изучается проблема уравнения Пуассона на квадратной области с использованием периодических граничных условий .Следовательно, температура вблизи границ не может уменьшаться до нуля.

Это причина, по которой нулевая частота, то есть среднее значение исходного члена, не используется.Если бы это было не так, стационарное уравнение теплопроводности не могло бы быть решено.Действительно, из-за периодических граничных условий тепловой поток, выходящий с одной стороны, проникает с противоположной стороны.Как и в случае, когда область полностью изолирована, температура будет повышаться и никогда не достигнет стационарного решения.

Чтобы применить нулевое граничное решение Дирихле, необходимо использовать одно из преобразований синуса вместоПреобразование Фурье.См. Например:

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

Вот ответ (мой!), Показывающий фрагмент кода, решающий двумерное уравнение Пуассона, нулевое граничное условие Дирихле, посредством синусоидального преобразования FFTW,Это c ++, но перевести его на c просто.

fftw3 для пуассона с граничным условием Дирихле для всех сторон вычислительной области

Чтобы восстановить что-то близкое к радиальной симметрииПоддержка источника должна быть далека от краев домена.Лучшим решением может быть переключение на одномерное радиальное уравнение, если вы знаете, что область, граничные условия и термин источника соответствуют радиальной симметрии.

Answer 2

Я думаю, что одна проблема здесь:

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

k=0 на (0,0), а не (N/2,N/2).Вот как определяется DFT.k работает от 0 до N-1. Вы вводите фазовый сдвиг в пространственной области.Я не знаю, как это повлияет на ваши последующие вычисления. Вы пытаетесь центрировать спектр, умножая входное значение на pow(-1,i).Этот трюк работает только для четных N!

Если лучше правильно вычислить ваши kx и ky:

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

Это будет работать длялюбой N.Также обратите внимание, что мое определение k выше соответствует радиальной частоте ( ω ), а не целому числу k ДПФ.Это также отличается от вашего, я думаю (если это так, ваш неверный!).

Другая проблема заключается в том, что DFT не FT.ДПФ вычисляется в ограниченной пространственной области, которая считается периодической.Этот домен является квадратом.Трудно делать вращательно-симметричные вещи на квадратной сетке.Частота Найквиста (максимально возможная частота) в 1030 раз выше по диагонали по сравнению с осями.И период шаблона также в 1031 раз выше по диагонали.Таким образом, вы получаете больше псевдонимов вдоль осей, и ваш шаблон с большей вероятностью будет перекрывать себя вдоль осей.

Если вышеприведенное вызывает эффект, который вы видите, попробуйте увеличить размер вашего пространственного домена,и / или плотность выборки.