Enforce reality for whole convolution term

2 files changed, 5 insertions, 8 deletions

diff --git a/src/main.c b/src/main.c
index 3d8867e..3c617a3 100644
--- a/src/main.c
+++ b/src/main.c
@@ -77,7 +77,7 @@ int main (
     energy(K1, K2, N1, N2, nsteps, nu, delta, L, g, print_freq, nthreads);
   }
   else if(command==COMMAND_ENSTROPHY){
-    energy(K1, K2, N1, N2, nsteps, nu, delta, L, g, print_freq, nthreads);
+    enstrophy(K1, K2, N1, N2, nsteps, nu, delta, L, g, print_freq, nthreads);
   }
   else if(command==COMMAND_QUIET){
     quiet(K1, K2, N1, N2, nsteps, nu, delta, L, g, nthreads);
diff --git a/src/navier-stokes.c b/src/navier-stokes.c
index e553022..fe2951f 100644
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -100,7 +100,6 @@ int energy(
     ins_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
     
     if(t%print_freq==0){
-
       energy=0.;
       for(kx=-K1;kx<=K1;kx++){
 	for (ky=-K2;ky<=K2;ky++){
@@ -334,7 +333,6 @@ int ns_init_u(
   return 0;
 }
 
-
 // next time step for Irreversible Navier-Stokes equation
 int ins_step(
   _Complex double* u,
@@ -451,7 +449,8 @@ int ins_rhs(
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
       if(kx!=0 || ky!=0){
-        out[klookup(kx,ky,2*K1+1,2*K2+1)]=-4*M_PI*M_PI/L/L*nu*(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]+(*g)(kx,ky)+4*M_PI*M_PI/L/L/sqrt(kx*kx+ky*ky)*ifft.fft[klookup(kx,ky,N1,N2)];
+	// enforce the reality of u by adding ifft.fft(k) and the conjugate of ifft.fft(-k)
+        out[klookup(kx,ky,2*K1+1,2*K2+1)]=-4*M_PI*M_PI/L/L*nu*(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]+(*g)(kx,ky)+4*M_PI*M_PI/L/L/sqrt(kx*kx+ky*ky)*(ifft.fft[klookup(kx,ky,N1,N2)]+conj(ifft.fft[klookup(-kx,-ky,N1,N2)]))/2;
       }
     }
   }
@@ -495,8 +494,7 @@ int ns_T(
   fftw_execute(fft2.fft_plan);
   // write to ifft
   for(i=0;i<N1*N2;i++){
-    // control numerical truncation by taking imaginary part (fft1 and fft2 should be purely imaginary)
-    ifft.fft[i]=-(__imag__ fft1.fft[i])*(__imag__ fft2.fft[i]);
+    ifft.fft[i]=fft1.fft[i]*fft2.fft[i];
   }
 
   // F(py/|p|*u)*F(qx*|q|*u)
@@ -520,8 +518,7 @@ int ns_T(
   fftw_execute(fft2.fft_plan);
   // write to ifft
   for(i=0;i<N1*N2;i++){
-    // control numerical truncation by taking imaginary part (fft1 and fft2 should be purely imaginary)
-    ifft.fft[i]=ifft.fft[i]+(__imag__ fft1.fft[i])*(__imag__ fft2.fft[i]);
+    ifft.fft[i]=ifft.fft[i]-fft1.fft[i]*fft2.fft[i];
   }
 
   // inverse fft