I am trying to numerically solve the Swift-Hohenberg equation http://en.wikipedia.org/wiki/Swift%E2%80%93Hohenberg_equation using a pseudo-spectral scheme, where the linear terms are treated implicitly in Fourier space, while the nonlinearity is evaluated in real space. A simple Euler scheme is used for the time integration.
My problem is that the Matlab code I have come up with works perfectly, while the C++ code, which relies on FFTW for the Fourier transforms, becomes unstable and diverges after a couple thousand time steps. I have tracked it down to the way the nonlinear term is treated (see the comments in the C++ code). If I use only the real part of Phi, the instability occurs. However, Phi should only have a negligible imaginary part due to numerical rounding errors, and Matlab is doing something similar, keeping Phi purely real.
The Matlab code also runs fine under Octave. The initial condition can be something like
R=0.02*(rand(256,256)-0.5);
in Matlab (small amplitude fluctuations).
TLDR;
Why do these pieces of code do different things? Specifically, how can I make the C++ code work the same way the Matlab version does?
Edit 1:
For completeness, I have added the code using the R2C/C2R functions provided by FFTW. See http://fftw.org/fftw3_doc/Multi_002dDimensional-DFTs-of-Real-Data.html for details (I hope I got the data layout right). This code always shows the instability after about 3100 time steps. If I reduce dt to e.g. 0.01, it occurs 10 times later.
C++ code using complex DFTs
#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++ code using 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 code
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
