Ich kämpfe seit fast einem Monat mit einem sehr seltsamen Fehler. Euch zu fragen ist meine letzte Hoffnung. Ich habe ein Programm in C geschrieben, das die 2d integriert Cahn-Hilliard-Gleichung unter Verwendung des impliziten Euler-Schemas (IE) im Fourier- (oder reziproken) Raum:
Wobei die “Hüte” bedeuten, dass wir uns im Fourier-Raum befinden: h_q(t_n+1) und h_q(t_n) sind die FTs von h(x,y) zu den Zeiten t_n und t_(n+1), N[h_q] ist der nichtlineare Operator, der auf h_q im Fourier-Raum angewendet wird, und L_q ist der lineare Operator, wiederum im Fourier-Raum. Ich möchte nicht zu sehr auf die Details der von mir verwendeten numerischen Methode eingehen, da ich sicher bin, dass das Problem nicht von dort kommt (ich habe versucht, andere Schemata zu verwenden).
Mein Code ist eigentlich ganz einfach. Hier ist der Anfang, wo ich im Wesentlichen Variablen deklariere, Speicher alloziiere und die Pläne für die FFTW-Routinen erstelle.
# 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);
Nach dieser Präambel initialisiere ich meine Funktion h(x,y) mit einem gleichmäßigen Zufallsrauschen und nehme auch die FT davon. Ich setze den Imaginärteil von h(x,y), also dh[i*Ny+j][1] im Code auf 0, da es sich um eine echte Funktion handelt. Dann berechne ich die Wellenvektoren qx und qyund mit ihnen berechne ich den linearen Operator meiner Gleichung im Fourier-Raum, der ist Linft im Code. Ich betrachte nur die -vierte Ableitung von h als linearen Term, sodass die FT des linearen Terms einfach -q^4 ist … aber ich möchte auch hier nicht auf die Details meiner Integrationsmethode eingehen. Die Frage bezieht sich nicht darauf.
// 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];
}
}
Dann kommt endlich die Zeitschleife. Im Wesentlichen mache ich folgendes:
-
Hin und wieder speichere ich die Daten in einer Datei und drucke einige Informationen auf dem Terminal aus. Insbesondere drucke ich den höchsten Wert der FT des nichtlinearen Terms. Ich überprüfe auch, ob h(x,y) gegen unendlich divergiert (das sollte nicht passieren!),
-
Berechnen Sie h^3 im direkten Raum (das ist einfach
dh[i*Ny+j][0]*dh[i*Ny+j][0]*dh[i*Ny+j][0]). Auch hier wird der Imaginärteil auf 0 gesetzt, -
Nimm die FT von h^3,
-
Ermitteln Sie den vollständigen nichtlinearen Term im reziproken Raum (d. h. N[h_q] im oben beschriebenen IE-Algorithmus) durch Berechnen von -q^2*(FT[h^3] -FT[h]). Im Code beziehe ich mich auf die Zeilen
Nonlft[i*Ny+j][0] = -Q2[i*Ny+j]*(Nonlft[i*Ny+j][0] -dhft[i*Ny+j][0])und die untere für den Imaginärteil. Ich mache das weil:
- Mit der IE-Methode in der Zeit vorrücken, in den direkten Raum zurücktransformieren und dann normalisieren.
Hier ist der Code:
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);
}
}
}
Letzter Teil des Codes: Speicher leeren und FFTW-Pläne zerstören.
// 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;
}
Wenn ich diesen Code ausführe, erhalte ich die folgende Ausgabe:
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!
Der Code stürzt ab, bevor er das Ende erreicht, und wir können sehen, dass der nichtlineare Term divergiert.
Was für mich keinen Sinn ergibt, ist, dass ich die Zeilen, in denen ich die FT des nichtlinearen Terms berechne, folgendermaßen ändere:
// 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];
}
}
Das bedeutet, dass ich diese Definition verwende:
statt dieser hier:
Dann ist der Code absolut stabil und es kommt zu keiner Divergenz! Sogar für Milliarden von Zeitschritten! Warum geschieht dies, da die zwei Arten der Berechnung Nonlft sollte gleichwertig sein?
Vielen Dank an alle, die sich die Zeit nehmen, das alles zu lesen und mir etwas zu helfen!
BEARBEITEN: Um die Dinge noch seltsamer zu machen, sollte ich darauf hinweisen, dass dieser Fehler NICHT für dasselbe System in 1D auftritt. In 1D beide Berechnungsmethoden Nonlft sind stabil.
BEARBEITEN: Ich füge eine kurze Animation hinzu, die zeigt, was mit der Funktion h(x,y) kurz vor dem Absturz passiert. Außerdem: Ich habe den Code in MATLAB schnell neu geschrieben, das Fast-Fourier-Transformationsfunktionen basierend auf der FFTW-Bibliothek verwendet, und der Fehler tritt NICHT auf … das Rätsel vertieft sich.
Kommentare sind nicht für längere Diskussionen gedacht; diese Konversation wurde in den Chat verschoben.
– Benutzer229044
♦
7. Dezember 2018 um 16:42 Uhr