/************************************************************************************************ * * * Ginz2D.c * * * * * * Solves the 2D Ginzburg-Landau-equation using a pseudospectral scheme in spatial domain * * and a exponential time-differencing scheme in temporal domain. * * The program displays its results using dislin. * * * * * * Some parameter sets: * * (alpha, beta) = ( 1.0, 2.0) //plain waves * * = ( 2.0, -1.0) //phase turbulence * * = ( 2.0, -2.0) //defect turbulence * * = ( 0.5, -1.5) //intermittency * * = ( 0.0, -4.0) //intermittency * * = ( 0.0, 1.5) //coherent structures, spirals * * * * These parameters are actually borrowed from the 1D simulation, but it shows * * that the parameters produce similar structures in 2D. * * * * * * compile with * * gcc -std=c99 -O3 -Wall -o Ginz2D Ginz2D.c -I(FFTW3_PATH)/include * * -L(FFTW_PATH)/lib -lfftw3 -I(DISLIN_PATH) -L(DISLIN_PATH) -ldislin -lm * * * * * * author: Johannes Luelff (johannes.luelff@uni-muenster.de) * * * ************************************************************************************************/ #include <stdio.h> #include <stdlib.h> //we want to use the C99 complex data type, so we need complex.h #include <complex.h> #include <math.h> #include "fftw3.h" #include "dislin.h" //these two defines control whether we... //... want to use dislin (we surely want to, cause otherwise // this program wouldnt produce any output) #define USE_DISLIN //... want dislin to print its result into a PNG file or onto the screen #undef USE_PNG //constants etc. #define PI 3.14159265358979323846264338327950288419716939937510582097494459230781 const int N = 256; //Grid resolution const double L = 200; //Physical length const double alpha = 0.0; //coefficients of GL-equation const double beta = 1.5; // -"- const int MAX_ITER = 1000; //maximal number of iterations; const int EVERY_DISPLAY = 5; //display every 'Every_DISPLAY'-th timestep //time(step) double t = 0; double dt = 0.05; //arrays etc. complex *A; //the 'main' field complex *Aedt_N1; //field for nonlinearity... complex *Aedt_N2; //... and backup field for nonlin. from previous timestep complex *Atmp; double *kabsvec; //vector that stores |k|**2 float *A_dislin_left, *A_dislin_right; //used for display with dislin fftw_plan planS2F, planF2S; //plan for spatial to fourier, fourier to spatial //multiplies array with scalar void smultArray(complex *arr, double s) { for(int i=0; i<N*N; i++) arr[i] *= s; } //fft from spatial to fourier domain void fftS2F(complex *arr) { fftw_execute_dft(planS2F, arr, arr); smultArray(arr, 1/sqrt(N*N)); } //fft from fourier to spatial domain void fftF2S(complex *arr) { fftw_execute_dft(planF2S, arr, arr); smultArray(arr, 1/sqrt(N*N)); } //assigns one array to the other void assignArrays(complex *out, const complex *in) { for(int i=0; i<N*N; i++) out[i] = in[i]; } //returns |z|**2 double cnorm(const complex z) { return creal(z)*creal(z) + cimag(z)*cimag(z); } //calculates the nonlinear part of the GL-equation in spatial domain void calcNonlin(const complex *in, complex *out) { assignArrays(out, in); fftF2S(out); for(int i=0; i<N*N; i++) out[i] = -(1.0+I*beta) * cnorm(out[i]) * out[i]; fftS2F(out); } //calculates the next time step using a 2nd order exponential time differencing scheme void doTimeStep_EDT2() { //switch nonlinearities N1 / N2 complex *ctmp = Aedt_N1; Aedt_N1 = Aedt_N2; Aedt_N2 = ctmp; //calc new nonlinearity calcNonlin(A, Aedt_N1); //put all together according to EDT2 //this may look scary, but its just the formula for EDT2 in a quite straight-forward manner. //q is the factor of the linear part of the GL-eq. //since we use a pseudospectral scheme, the 2nd derivative is also just a factor for(int i=0; i<N*N; i++) { complex q = (-(1+I*alpha)*kabsvec[i]+1); A[i] = A[i] * cexp(q*dt) + Aedt_N1[i] * ((1+dt*q)*cexp(dt*q) - 1 - 2*dt*q) / (dt*q*q) + Aedt_N2[i] * (-cexp(dt*q)+1+dt*q) / (dt*q*q); } t += dt; } //displays real part and absolute value of the field using dislin //dont take all the numbers regarding dislin in this routine as 'the only truth'; //dealing with dislin is always kind of tricky and try-and-error-like void display(complex *in) { //to spatial domain assignArrays(Atmp, in); fftF2S(Atmp); //prepare output array for dislin (left panel) double min_left = HUGE_VAL; double max_left = -HUGE_VAL; for(int i=0; i<N*N; i++) { A_dislin_left[i] = creal(Atmp[i]); if(A_dislin_left[i]<min_left) min_left = A_dislin_left[i]; if(A_dislin_left[i]>max_left) max_left = A_dislin_left[i]; } if(fabs(min_left-max_left)<1e-6) { min_left=-0.1; max_left= 0.1; } //prepare output array for dislin (right panel) double min_right = HUGE_VAL; double max_right = -HUGE_VAL; for(int i=0; i<N*N; i++) { A_dislin_right[i] = cabs(Atmp[i]); if(A_dislin_right[i]<min_right) min_right = A_dislin_right[i]; if(A_dislin_right[i]>max_right) max_right = A_dislin_right[i]; } if(fabs(min_right-max_right)<1e-6) { min_right=-0.1; max_right= 0.1; } //plot with dislin #ifdef USE_PNG //init dislin metafl("PNG"); setfil("out.png"); scrmod("REVERS"); winsiz(1280,640); page(2970*2,2970); disini(); ax3len(2970*0.75, 2970*0.75, 2970*0.75); #endif //USE_PNG erase(); axspos(2970*0.1,2970*0.85); graf3(0,L,0,L/2, 0,L,0,L/2, min_left, max_left, min_left, (max_left-min_left)/2); crvmat(A_dislin_left, N, N, ((512/N)<=0)?1:(512/N), ((512/N)<=0)?1:(512/N)); titlin("Real(A)", 3); title(); endgrf(); axspos(2970*1.1,2970*0.85); graf3(0,L,0,L/2, 0,L,0,L/2, min_right, max_right, min_right, (max_right-min_right)/2); crvmat(A_dislin_right, N, N, ((512/N)<=0)?1:(512/N), ((512/N)<=0)?1:(512/N)); titlin("Abs(A)", 3); title(); endgrf(); #ifdef USE_PNG disfin(); #endif //USE_PNG } int main(int argc, char *argv[]) { //get memory A = fftw_malloc(sizeof(complex)*N*N); Atmp = fftw_malloc(sizeof(complex)*N*N); Aedt_N1 = fftw_malloc(sizeof(complex)*N*N); Aedt_N2 = fftw_malloc(sizeof(complex)*N*N); kabsvec = fftw_malloc(sizeof(double)*N*N); A_dislin_left = fftw_malloc(sizeof(float)*N*N); A_dislin_right = fftw_malloc(sizeof(float)*N*N); //make plans for fftw planS2F = fftw_plan_dft_2d(N, N, (fftw_complex*)A, (fftw_complex*)A, FFTW_FORWARD, FFTW_ESTIMATE); planF2S = fftw_plan_dft_2d(N, N, (fftw_complex*)A, (fftw_complex*)A, FFTW_BACKWARD, FFTW_ESTIMATE); //init kvec for(int i=0; i<N; i++) for(int j=0; j<N; j++) { double kx = 2*PI * ((i<N/2+1) ? (i) : (i-N)) / L; double ky = 2*PI * ((j<N/2+1) ? (j) : (j-N)) / L; kabsvec[i*N+j] = kx*kx + ky*ky; } //init field with white noise of amplitude 0.01 //and go to fourier space, cause we will do the timestepping in fourier space for(int i=0; i<N*N; i++) A[i] = 0.02*(0.5 - rand()/(double)RAND_MAX + I*(0.5 - rand()/(double)RAND_MAX)); fftS2F(A); #ifdef USE_DISLIN #ifndef USE_PNG //init dislin metafl("XWIN"); scrmod("REVERS"); winsiz(1280,640); page(2970*2,2970); disini(); ax3len(2970*0.75, 2970*0.75, 2970*0.75); #endif //USE_PNG #endif //USE_DISLIN //main timestepping loop for(int iter=0; iter<=MAX_ITER; iter++) { doTimeStep_EDT2(); #ifdef USE_DISLIN if(iter%EVERY_DISPLAY==0) { display(A); printf("t = %e , dt = %e\n", t, dt); } #endif //USE_DISLIN } #ifdef USE_DISLIN #ifndef USE_PNG //finalize dislin disfin(); #endif //USE_PNG #endif //USE_DISLIN //clean up stuff free(A); free(Atmp); free(Aedt_N1); free(Aedt_N2); free(kabsvec); free(A_dislin_left); free(A_dislin_right); fftw_destroy_plan(planF2S); fftw_destroy_plan(planS2F); //yay were done printf("READY\n"); }