Почти месяц бьюсь с очень странным багом. Спросить вас, ребята, моя последняя надежда. Я написал программу на C, которая интегрирует двумерное уравнение Кана–Хиллиарда, используя Неявная схема Эйлера (IE) в пространстве Фурье (или обратном):
Где «шляпы» означают, что мы находимся в пространстве Фурье: h_q(t_n+1) и h_q(t_n) — FTs h(x,y) в моменты времени t_n и t_(n+1), N[h_q] — нелинейный оператор применяется к h_q в пространстве Фурье, а L_q является линейным, опять же в пространстве Фурье. Я не хочу слишком вдаваться в подробности используемого численного метода, так как уверен, что проблема не оттуда (пробовал использовать другие схемы).
Мой код на самом деле довольно прост. Вот начало, где в основном я объявляю переменные, выделяю память и создаю планы для подпрограмм FFTW.
# include <stdlib.h>
# include <stdio.h>
# include <time.h>
# include <math.h>
# include <fftw3.h>
# define pi M_PI
int main(){
// define lattice size and spacing
int Nx = 150; // n of points on x
int Ny = 150; // n of points on y
double dx = 0.5; // bin size on x and y
// define simulation time and time step
long int Nt = 1000; // n of time steps
double dt = 0.5; // time step size
// number of frames to plot (at denominator)
long int nframes = Nt/100;
// define the noise
double rn, drift = 0.05; // punctual drift of h(x)
srand(666); // seed the RNG
// other variables
int i, j, nt; // variables for space and time loops
// declare FFTW3 routine
fftw_plan FT_h_hft; // routine to perform fourier transform
fftw_plan FT_Nonl_Nonlft;
fftw_plan IFT_hft_h; // routine to perform inverse fourier transform
// declare and allocate memory for real variables
double *Linft = fftw_alloc_real(Nx*Ny);
double *Q2 = fftw_alloc_real(Nx*Ny);
double *qx = fftw_alloc_real(Nx);
double *qy = fftw_alloc_real(Ny);
// declare and allocate memory for complex variables
fftw_complex *dh = fftw_alloc_complex(Nx*Ny);
fftw_complex *dhft = fftw_alloc_complex(Nx*Ny);
fftw_complex *Nonl = fftw_alloc_complex(Nx*Ny);
fftw_complex *Nonlft = fftw_alloc_complex(Nx*Ny);
// create the FFTW plans
FT_h_hft = fftw_plan_dft_2d ( Nx, Ny, dh, dhft, FFTW_FORWARD, FFTW_ESTIMATE );
FT_Nonl_Nonlft = fftw_plan_dft_2d ( Nx, Ny, Nonl, Nonlft, FFTW_FORWARD, FFTW_ESTIMATE );
IFT_hft_h = fftw_plan_dft_2d ( Nx, Ny, dhft, dh, FFTW_BACKWARD, FFTW_ESTIMATE );
// open file to store the data
char acstr[160];
FILE *fp;
sprintf(acstr, "CH2d_IE_dt%.2f_dx%.3f_Nt%ld_Nx%d_Ny%d_#f%.ld.dat",dt,dx,Nt,Nx,Ny,Nt/nframes);
После этой преамбулы я инициализирую свою функцию h(x,y) равномерным случайным шумом, а также беру для нее Фурье-Фон. Я установил мнимую часть h(x,y), которая в коде равна dh[i*Ny+j][1], равной 0, так как это реальная функция. Затем я вычисляю волновые векторы qx и qy и с их помощью вычисляю линейный оператор моего уравнения в пространстве Фурье, который в коде равен Linft. Я рассматриваю только - четвертую производную от h как линейный член, так что FT линейного члена просто -q ^ 4... но опять же, я не хочу вдаваться в подробности моего метода интегрирования. Вопрос не в этом.
// generate h(x,y) at initial time
for ( i = 0; i < Nx; i++ ) {
for ( j = 0; j < Ny; j++ ) {
rn = (double) rand()/RAND_MAX; // extract a random number between 0 and 1
dh[i*Ny+j][0] = drift-2.0*drift*rn; // shift of +-drift
dh[i*Ny+j][1] = 0.0;
}
}
// execute plan for the first time
fftw_execute (FT_h_hft);
// calculate wavenumbers
for (i = 0; i < Nx; i++) { qx[i] = 2.0*i*pi/(Nx*dx); }
for (i = 0; i < Ny; i++) { qy[i] = 2.0*i*pi/(Ny*dx); }
for (i = 1; i < Nx/2; i++) { qx[Nx-i] = -qx[i]; }
for (i = 1; i < Ny/2; i++) { qy[Ny-i] = -qy[i]; }
// calculate the FT of the linear operator
for ( i = 0; i < Nx; i++ ) {
for ( j = 0; j < Ny; j++ ) {
Q2[i*Ny+j] = qx[i]*qx[i] + qy[j]*qy[j];
Linft[i*Ny+j] = -Q2[i*Ny+j]*Q2[i*Ny+j];
}
}
Затем, наконец, наступает петля времени. По сути, я делаю следующее:
Время от времени я сохраняю данные в файл и распечатываю некоторую информацию на терминале. В частности, я печатаю наибольшее значение FT нелинейного термина. Я также проверяю, расходится ли h(x,y) в бесконечность (этого не должно быть!),
Вычислите h^3 в прямом пространстве (то есть просто
dh[i*Ny+j][0]*dh[i*Ny+j][0]*dh[i*Ny+j][0]). Опять же, мнимая часть установлена на 0,Возьмите FT h ^ 3,
Получите полный нелинейный член в обратном пространстве (то есть N[h_q] в алгоритме IE, написанном выше), вычислив -q^2*(FT[h^3] - FT[h]). В коде я имею в виду строки
Nonlft[i*Ny+j][0] = -Q2[i*Ny+j]*(Nonlft[i*Ny+j][0] -dhft[i*Ny+j][0])и одну ниже для мнимой части. Я делаю это, потому что:
- Продвиньтесь во времени, используя метод IE, преобразуйте обратно в прямое пространство, а затем нормализуйте.
Вот код:
for(nt = 0; nt < Nt; nt++) {
if((nt % nframes)== 0) {
printf("%.0f %%\n",((double)nt/(double)Nt)*100);
printf("Nonlft %.15f \n",Nonlft[(Nx/2)*(Ny/2)][0]);
// write data to file
fp = fopen(acstr,"a");
for ( i = 0; i < Nx; i++ ) {
for ( j = 0; j < Ny; j++ ) {
fprintf(fp, "%4d %4d %.6f\n", i, j, dh[i*Ny+j][0]);
}
}
fclose(fp);
}
// check if h is going to infinity
if (isnan(dh[1][0])!=0) {
printf("crashed!\n");
return 0;
}
// calculate nonlinear term h^3 in direct space
for ( i = 0; i < Nx; i++ ) {
for ( j = 0; j < Ny; j++ ) {
Nonl[i*Ny+j][0] = dh[i*Ny+j][0]*dh[i*Ny+j][0]*dh[i*Ny+j][0];
Nonl[i*Ny+j][1] = 0.0;
}
}
// Fourier transform of nonlinear term
fftw_execute (FT_Nonl_Nonlft);
// second derivative in Fourier space is just multiplication by -q^2
for ( i = 0; i < Nx; i++ ) {
for ( j = 0; j < Ny; j++ ) {
Nonlft[i*Ny+j][0] = -Q2[i*Ny+j]*(Nonlft[i*Ny+j][0] -dhft[i*Ny+j][0]);
Nonlft[i*Ny+j][1] = -Q2[i*Ny+j]*(Nonlft[i*Ny+j][1] -dhft[i*Ny+j][1]);
}
}
// Implicit Euler scheme in Fourier space
for ( i = 0; i < Nx; i++ ) {
for ( j = 0; j < Ny; j++ ) {
dhft[i*Ny+j][0] = (dhft[i*Ny+j][0] + dt*Nonlft[i*Ny+j][0])/(1.0 - dt*Linft[i*Ny+j]);
dhft[i*Ny+j][1] = (dhft[i*Ny+j][1] + dt*Nonlft[i*Ny+j][1])/(1.0 - dt*Linft[i*Ny+j]);
}
}
// transform h back in direct space
fftw_execute (IFT_hft_h);
// normalize
for ( i = 0; i < Nx; i++ ) {
for ( j = 0; j < Ny; j++ ) {
dh[i*Ny+j][0] = dh[i*Ny+j][0] / (double) (Nx*Ny);
dh[i*Ny+j][1] = dh[i*Ny+j][1] / (double) (Nx*Ny);
}
}
}
Последняя часть кода: очистить память и уничтожить планы FFTW.
// terminate the FFTW3 plan and free memory
fftw_destroy_plan (FT_h_hft);
fftw_destroy_plan (FT_Nonl_Nonlft);
fftw_destroy_plan (IFT_hft_h);
fftw_cleanup();
fftw_free(dh);
fftw_free(Nonl);
fftw_free(qx);
fftw_free(qy);
fftw_free(Q2);
fftw_free(Linft);
fftw_free(dhft);
fftw_free(Nonlft);
return 0;
}
Если я запускаю этот код, я получаю следующий вывод:
0 %
Nonlft 0.0000000000000000000
1 %
Nonlft -0.0000000000001353512
2 %
Nonlft -0.0000000000000115539
3 %
Nonlft 0.0000000001376379599
...
69 %
Nonlft -12.1987455309071730625
70 %
Nonlft -70.1631962517720353389
71 %
Nonlft -252.4941743351609204637
72 %
Nonlft 347.5067875825179726235
73 %
Nonlft 109.3351142318568633982
74 %
Nonlft 39933.1054502610786585137
crashed!
Код вылетает, не дойдя до конца, и мы видим, что нелинейный член расходится.
Теперь то, что для меня не имеет смысла, это то, что если я изменю строки, в которых я вычисляю FT нелинейного термина, следующим образом:
// calculate nonlinear term h^3 -h in direct space
for ( i = 0; i < Nx; i++ ) {
for ( j = 0; j < Ny; j++ ) {
Nonl[i*Ny+j][0] = dh[i*Ny+j][0]*dh[i*Ny+j][0]*dh[i*Ny+j][0] -dh[i*Ny+j][0];
Nonl[i*Ny+j][1] = 0.0;
}
}
// Fourier transform of nonlinear term
fftw_execute (FT_Nonl_Nonlft);
// second derivative in Fourier space is just multiplication by -q^2
for ( i = 0; i < Nx; i++ ) {
for ( j = 0; j < Ny; j++ ) {
Nonlft[i*Ny+j][0] = -Q2[i*Ny+j]* Nonlft[i*Ny+j][0];
Nonlft[i*Ny+j][1] = -Q2[i*Ny+j]* Nonlft[i*Ny+j][1];
}
}
Это означает, что я использую это определение:
вместо этого:
Тогда код абсолютно стабилен и никаких расхождений не происходит! Даже за миллиарды временных шагов! Почему это происходит, если два способа расчета Nonlft должны быть эквивалентны?
Большое спасибо всем, кто найдет время, чтобы прочитать все это и помочь мне!
РЕДАКТИРОВАТЬ: Чтобы сделать вещи еще более странными, я должен указать, что эта ошибка НЕ происходит для той же системы в 1D. В 1D оба метода расчета Nonlft стабильны.
РЕДАКТИРОВАТЬ: я добавляю короткую анимацию того, что происходит с функцией h(x,y) непосредственно перед сбоем. Кроме того: я быстро переписал код в MATLAB, который использует функции быстрого преобразования Фурье, основанные на библиотеке FFTW, и ошибки НЕ происходит... тайна углубляется. 
