Simpler expression for fft term

1 files changed, 16 insertions, 22 deletions

diff --git a/src/navier-stokes.c b/src/navier-stokes.c
index 3dd7484..86a8497 100644
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -68,7 +68,7 @@ int ins_step(_Complex double* u, ns_params params, fft_vects vects, _Complex dou
 int ins_rhs(_Complex double* out, _Complex double* u, ns_params params, fft_vects vects){
   int kx,ky;
 
-  // F(u/|p|)*F(q1*q2*u/|q|)
+  // F(px/|p|*u)*F(qy*|q|*u)
   // init to 0
   for(kx=0; kx<params.N*params.N; kx++){
     vects.fft1[kx]=0;
@@ -79,11 +79,12 @@ int ins_rhs(_Complex double* out, _Complex double* u, ns_params params, fft_vect
   for(kx=-params.K;kx<=params.K;kx++){
     for(ky=-params.K;ky<=params.K;ky++){
       if(kx!=0 || ky!=0){
-        vects.fft1[KLOOKUP(kx,ky,params.N)]=u[KLOOKUP(kx,ky,params.S)]/sqrt(kx*kx+ky*ky);
-        vects.fft2[KLOOKUP(kx,ky,params.N)]=kx*ky*u[KLOOKUP(kx,ky,params.S)]/sqrt(kx*kx+ky*ky);
+        vects.fft1[KLOOKUP(kx,ky,params.N)]=kx/sqrt(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)];
+        vects.fft2[KLOOKUP(kx,ky,params.N)]=ky*sqrt(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)];
       }
     }
   }
+
   // fft
   fftw_execute(vects.fft1_plan);
   fftw_execute(vects.fft2_plan);
@@ -93,51 +94,44 @@ int ins_rhs(_Complex double* out, _Complex double* u, ns_params params, fft_vect
       vects.invfft[KLOOKUP(kx,ky,params.N)]=vects.fft1[KLOOKUP(kx,ky,params.N)]*vects.fft2[KLOOKUP(kx,ky,params.N)];
     }
   }
-  // inverse fft
-  fftw_execute(vects.invfft_plan);
-  
-  // write out
-  for(kx=0; kx<params.S*params.S; kx++){
-    out[kx]=0;
-  }
-  for(kx=-params.K;kx<=params.K;kx++){
-    for(ky=-params.K;ky<=params.K;ky++){
-      if(kx!=0 || ky!=0){
-        out[KLOOKUP(kx,ky,params.S)]=-4*M_PI*M_PI*params.nu*(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)]+params.g[KLOOKUP(kx,ky,params.S)]+4*M_PI*M_PI/sqrt(kx*kx+ky*ky)*vects.invfft[KLOOKUP(kx,ky,params.N)]*(kx*kx-ky*ky)/params.N/params.N;
-      }
-    }
-  }
 
-  // F(u/|p|)*F((q1*q1-q2*q2)*u/|q|)
+  // F(py/|p|*u)*F(qx*|q|*u)
   // init to 0
   for(kx=0; kx<params.N*params.N; kx++){
+    vects.fft1[kx]=0;
     vects.fft2[kx]=0;
-    vects.invfft[kx]=0;
   }
   // fill modes
   for(kx=-params.K;kx<=params.K;kx++){
     for(ky=-params.K;ky<=params.K;ky++){
       if(kx!=0 || ky!=0){
-        vects.fft2[KLOOKUP(kx,ky,params.N)]=(kx*kx-ky*ky)*u[KLOOKUP(kx,ky,params.S)]/sqrt(kx*kx+ky*ky);
+        vects.fft1[KLOOKUP(kx,ky,params.N)]=ky/sqrt(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)];
+        vects.fft2[KLOOKUP(kx,ky,params.N)]=kx*sqrt(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)];
       }
     }
   }
+
   // fft
+  fftw_execute(vects.fft1_plan);
   fftw_execute(vects.fft2_plan);
   // write to invfft
   for(kx=-2*params.K;kx<=2*params.K;kx++){
     for(ky=-2*params.K;ky<=2*params.K;ky++){
-      vects.invfft[KLOOKUP(kx,ky,params.N)]=vects.fft1[KLOOKUP(kx,ky,params.N)]*vects.fft2[KLOOKUP(kx,ky,params.N)];
+      vects.invfft[KLOOKUP(kx,ky,params.N)]=vects.invfft[KLOOKUP(kx,ky,params.N)]-vects.fft1[KLOOKUP(kx,ky,params.N)]*vects.fft2[KLOOKUP(kx,ky,params.N)];
     }
   }
+
   // inverse fft
   fftw_execute(vects.invfft_plan);
   
   // write out
+  for(kx=0; kx<params.S*params.S; kx++){
+    out[kx]=0;
+  }
   for(kx=-params.K;kx<=params.K;kx++){
     for(ky=-params.K;ky<=params.K;ky++){
       if(kx!=0 || ky!=0){
-        out[KLOOKUP(kx,ky,params.S)]=out[KLOOKUP(kx,ky,params.S)]-4*M_PI*M_PI/sqrt(kx*kx+ky*ky)*vects.invfft[KLOOKUP(kx,ky,params.N)]*(kx*ky)/params.N/params.N;
+        out[KLOOKUP(kx,ky,params.S)]=-4*M_PI*M_PI*params.nu*(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)]+params.g[KLOOKUP(kx,ky,params.S)]+4*M_PI*M_PI/sqrt(kx*kx+ky*ky)*vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N;
       }
     }
   }