Численная нестабильность FFTW <> Matlab - PullRequest
Новая
       23

Численная нестабильность FFTW <> Matlab

9 голосов
/ 27 июля 2011

Я пытаюсь численно решить уравнение Свифта-Хоенберга http://en.wikipedia.org/wiki/Swift%E2%80%93Hohenberg_equation, используя псевдоспектральную схему, где линейные члены обрабатываются неявно в пространстве Фурье, а нелинейность оценивается в реальном пространстве.Простая схема Эйлера используется для временной интеграции.
Моя проблема в том, что код Matlab, который я придумал, работает отлично, в то время как код C ++, который использует FFTW для преобразований Фурье, становится нестабильным и расходится после парытысяча временных шагов.Я проследил это до того, как трактуется нелинейный термин (см. Комментарии в коде C ++).Если я использую только настоящую часть Фи, возникает нестабильность.Однако из-за ошибок округления чисел Phi должна иметь лишь незначительную мнимую часть, а Matlab делает нечто подобное, сохраняя Phi чисто реальным.Код Matlab также отлично работает под Octave.Начальное условие может быть что-то вроде
R=0.02*(rand(256,256)-0.5);
в Matlab (небольшие колебания амплитуды).

TLDR;

Почему эти части кода делают разные вещи?В частности, как я могу сделать так, чтобы код C ++ работал так же, как и версия Matlab?

Редактировать 1:

Для полноты я добавил код, используя функции R2C / C2R, предоставляемые FFTW.,См. http://fftw.org/fftw3_doc/Multi_002dDimensional-DFTs-of-Real-Data.html для деталей (надеюсь, я правильно понял расположение данных).Этот код всегда показывает нестабильность после примерно 3100 временных шагов.Если я уменьшу dt, например, до 0,01, это произойдет в 10 раз позже.

Код C ++ с использованием сложных кодов DFT 

#include <iostream>
#include <fstream>
#include <cmath>
#include <fftw3.h>

int main() {

const int N=256, nSteps=10000;
const double k=2.0*M_PI/N, dt=0.1, eps=0.25;

double *Buf=(double*)fftw_malloc(N*N*sizeof(double));
double *D0=(double*)fftw_malloc(N*N*sizeof(double));

// complex arrays
fftw_complex *Phi=(fftw_complex*)fftw_malloc(N*N*sizeof(fftw_complex));
fftw_complex *PhiF=(fftw_complex*)fftw_malloc(N*N*sizeof(fftw_complex));
fftw_complex *NPhiF=(fftw_complex*)fftw_malloc(N*N*sizeof(fftw_complex));

// plans for Fourier transforms
fftw_plan phiPlan=fftw_plan_dft_2d(N, N, Phi, PhiF, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_plan nPhiPlan=fftw_plan_dft_2d(N, N, NPhiF, NPhiF, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_plan phiInvPlan=fftw_plan_dft_2d(N, N, Phi, Phi, FFTW_BACKWARD, FFTW_ESTIMATE);

std::ifstream fin("R.dat", std::ios::in | std::ios::binary); // read initial condition
fin.read(reinterpret_cast<char*>(Buf), N*N*sizeof(double));
fin.close();
for(int i=0; i<N*N; i++) {
    Phi[i][0]=Buf[i];   //initial condition
    Phi[i][1]=0.0;  //no imaginary part
}

fftw_execute(phiPlan);  //PhiF contains FT of initial condition

for(int j=0; j<N; j++) {
    for(int i=0; i<N; i++) {

        double kx=(i-(i/(N-N/2)*N))*k;
        double ky=(j-(j/(N-N/2)*N))*k;
        double k2=kx*kx+ky*ky;

        D0[j*N+i]=1.0/((1.0 - dt*(eps-1.0 + 2.0*k2 - k2*k2)));  // array of prefactors
    }
}   

const double norm=1.0/(N*N);

for(int n=0; n<=nSteps; n++) {

    if(n%100==0) {
        std::cout<<"n = "<<n<<'\n';
    }

    for(int j=0; j<N*N; j++) {
        // nonlinear term Phi^3
        //NPhiF[j][0]=Phi[j][0]*Phi[j][0]*Phi[j][0]; // unstable
        //NPhiF[j][1]=0.0;
            NPhiF[j][0]=Phi[j][0]*Phi[j][0]*Phi[j][0] - 3.0*Phi[j][0]*Phi[j][1]*Phi[j][1];
            NPhiF[j][1]=-Phi[j][1]*Phi[j][1]*Phi[j][1] + 3.0*Phi[j][0]*Phi[j][0]*Phi[j][1];
    }

    fftw_execute(nPhiPlan); // NPhiF contains FT of Phi^3

    for(int j=0; j<N*N; j++) {
        PhiF[j][0]=(PhiF[j][0] - dt*NPhiF[j][0])*D0[j]; // update
        PhiF[j][1]=(PhiF[j][1] - dt*NPhiF[j][1])*D0[j];

        Phi[j][0]=PhiF[j][0]*norm; // FFTW does not normalize
        Phi[j][1]=PhiF[j][1]*norm;
    }

    fftw_execute(phiInvPlan); // Phi contains the updated Phi in real space
}

for(int i=0; i<N*N; i++) {
    Buf[i]=Phi[i][0];   // saving the real part of Phi
}
std::ofstream fout("Phi.dat", std::ios::trunc | std::ios::binary);
fout.write(reinterpret_cast<char*>(Buf), N*N*sizeof(double));
fout.close();

for(int i=0; i<N*N; i++) {
    Buf[i]=Phi[i][1];   // saving the imag part of Phi
}
fout.open("PhiImag.dat", std::ios::trunc | std::ios::binary);
fout.write(reinterpret_cast<char*>(Buf), N*N*sizeof(double));
fout.close();


fftw_free(D0);
fftw_free(Buf);

fftw_free(Phi);
fftw_free(PhiF);
fftw_free(NPhiF);

fftw_destroy_plan(phiPlan);
fftw_destroy_plan(phiInvPlan);
fftw_destroy_plan(nPhiPlan);

return EXIT_SUCCESS;
}

C ++ с использованием R2C / C2R 

#include <iostream>
#include <fstream>
#include <cmath>
#include <fftw3.h>

int main() {

const int N=256, nSteps=3100;
const int w=N/2+1;
const double k=2.0*M_PI/N, dt=0.1, eps=0.25;

double *Buf=(double*)fftw_malloc(N*N*sizeof(double));
double *D0=(double*)fftw_malloc(N*w*sizeof(double));

fftw_complex *Phi=(fftw_complex*)fftw_malloc(N*w*sizeof(fftw_complex));
fftw_complex *PhiF=(fftw_complex*)fftw_malloc(N*w*sizeof(fftw_complex));
fftw_complex *NPhi=(fftw_complex*)fftw_malloc(N*w*sizeof(fftw_complex));

fftw_plan phiPlan=fftw_plan_dft_r2c_2d(N, N, (double*)PhiF, PhiF, FFTW_ESTIMATE);
fftw_plan nPhiPlan=fftw_plan_dft_r2c_2d(N, N, (double*)NPhi, NPhi, FFTW_ESTIMATE);
fftw_plan phiInvPlan=fftw_plan_dft_c2r_2d(N, N, Phi, (double*)Phi, FFTW_ESTIMATE);

std::ifstream fin("R.dat", std::ios::in | std::ios::binary);
fin.read(reinterpret_cast<char*>(Buf), N*N*sizeof(double));
fin.close();
for(int j=0; j<N; j++) {
    for(int i=0; i<N; i++) {
        ((double*)PhiF)[j*2*w+i]=Buf[j*N+i];
        ((double*)Phi)[j*2*w+i]=Buf[j*N+i];
    }
}

fftw_execute(phiPlan); //PhiF contains FT of IC

for(int j=0; j<N; j++) {
    for(int i=0; i<w; i++) {

        double kx=(i-(i/(N-N/2)*N))*k;
        double ky=(j-(j/(N-N/2)*N))*k;
        double k2=kx*kx+ky*ky;

        D0[j*w+i]=1.0/(1.0 - dt*(eps-1.0 + 2.0*k2 - k2*k2));
    }
}

const double norm=1.0/(N*N);

//begin first Euler step
for(int n=0; n<=nSteps; n++) {

    if(n%100==0) {
        std::cout<<"n = "<<n<<'\n';
    }

    for(int j=0; j<N; j++) {
        for(int i=0; i<N; i++) {
            ((double*)NPhi)[j*2*w+i]=((double*)Phi)[j*2*w+i]  *((double*)Phi)[j*2*w+i]  * ((double*)Phi)[j*2*w+i];
        }
    }

    fftw_execute(nPhiPlan); // NPhi contains FT of Phi^3

    for(int j=0; j<N*w; j++) {
        PhiF[j][0]=(PhiF[j][0] - dt*NPhi[j][0])*D0[j];
        PhiF[j][1]=(PhiF[j][1] - dt*NPhi[j][1])*D0[j];
    }

    for(int j=0; j<N*w; j++) {
        Phi[j][0]=PhiF[j][0]*norm;
        Phi[j][1]=PhiF[j][1]*norm;
    }

    fftw_execute(phiInvPlan);

}

for(int j=0; j<N; j++) {
    for(int i=0; i<N; i++) {
        Buf[j*N+i]=((double*)Phi)[j*2*w+i];
    }
}

std::ofstream fout("Phi.dat", std::ios::trunc | std::ios::binary);
fout.write(reinterpret_cast<char*>(Buf), N*N*sizeof(double));
fout.close();

fftw_destroy_plan(phiPlan);
fftw_destroy_plan(phiInvPlan);
fftw_destroy_plan(nPhiPlan);

fftw_free(D0);
fftw_free(Buf);
fftw_free(Phi);
fftw_free(PhiF);
fftw_free(NPhi);
}

код Matlab 

function Phi=SwiHoEuler(Phi, nSteps)
epsi=0.25;
dt=0.1;

[nR nC]=size(Phi);
if mod(nR, 2)==0
    kR=[0:nR/2-1 -nR/2:-1]*2*pi/nR;
else
    kR=[0:nR/2 -floor(nR/2):-1]*2*pi/nR;
end
Ky=repmat(kR.', 1, nC);

if mod(nC, 2)==0
    kC=[0:nC/2-1 -nC/2:-1]*2*pi/nC;
else
    kC=[0:nC/2 -floor(nC/2):-1]*2*pi/nC;
end
Kx=repmat(kC, nR, 1); % frequencies
K2=Kx.^2+Ky.^2; % used for Laplacian in Fourier space
D0=1.0./(1.0-dt*(epsi-1.0+2.0*K2-K2.*K2)); % linear factors combined

PhiF=fft2(Phi);

for n=0:nSteps
    NPhiF=fft2(Phi.^3); % nonlinear term, evaluated in real space
    if mod(n, 100)==0
        fprintf('n = %i\n', n);
    end
    PhiF=(PhiF - dt*NPhiF).*D0; % update

    Phi=ifft2(PhiF); % inverse transform
end
return

Ответы [ 2 ]

1 голос
/ 23 января 2012

Посмотрите на эти строки: 

for ...
  double kx=(i-(i/(N-N/2)*N))*k;
  double ky=(j-(j/(N-N/2)*N))*k;
  double k2=kx*kx+ky*ky;
...

Я должен признать, что я не изучал алгоритмы, но "i / (N-N / 2)" состоит из целых чисел, и я подозреваю, что ваши kx, ky и k2 рассчитываются, как ожидалось. Вы можете попробовать что-то вроде этого, чтобы избежать потенциальных ошибок целочисленного округления: 

for ...
  double kx=( double(i) -( double(i)/(0.5*double(N*N)))*k; // where in our case: (N-N/2)*N) = 0.5*N*N
  ...
...
0 голосов
/ 27 июля 2011

РЕДАКТИРОВАТЬ Ниже приведен неправильный OP был прав.

Указатель на действительную часть хранится в [0], мнимый - в [1] (т. е. NPhi [1] [j] - это то, на что вы должны ссылаться - по крайней мере, согласно их сайту ). Так что это, вероятно, одна проблема. Фиксированные линии (я считаю) должны быть следующими: 

NPhiF[0][j]=Phi[0][j]*Phi[0][j]*Phi[0][j] - 3.0*Phi[0][j]*Phi[1][j]*Phi[1][j];
NPhiF[1][j]=-Phi[1][j]*Phi[1][j]*Phi[1][j] + 3.0*Phi[0][j]*Phi[0][j]*Phi[1][j];

Это исправит одну часть - вам придется исправить все остальное.

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.

Похожие темы

...