Simpler expression for fft term

2 files changed, 22 insertions, 26 deletions

diff --git a/src/main.c b/src/main.c
index 9420e90..b92698a 100644
--- a/src/main.c
+++ b/src/main.c
@@ -77,10 +77,10 @@ int read_args(int argc, const char* argv[], ns_params* params, unsigned int* nst
   *nsteps=16777216;
   params->nu=4.9632717887631524e-05;
   */
-  params->h=1e-5;
   params->K=16;
+  params->h=1e-3/(2*params->K+1);
   *nsteps=10000000;
-  params->nu=1e-4;
+  params->nu=1./1024/(2*params->K+1);
 
   // loop over arguments
   for(i=1;i<argc;i++){
@@ -221,7 +221,7 @@ int enstrophy(ns_params params, unsigned int Nsteps){
     for(ky=-params.K;ky<=params.K;ky++){
       //params.g[KLOOKUP(kx,ky,params.S)]=sqrt(kx*kx*ky*ky)*exp(-(kx*kx+ky*ky));
       if((kx==2 && ky==-1) || (kx==-2 && ky==1)){
-	params.g[KLOOKUP(kx,ky,params.S)]=1;
+	params.g[KLOOKUP(kx,ky,params.S)]=1.0*params.K;
       }
       else{
 	params.g[KLOOKUP(kx,ky,params.S)]=0;
@@ -246,19 +246,21 @@ int enstrophy(ns_params params, unsigned int Nsteps){
     ins_step(u, params, fft_vects, tmp1, tmp2, tmp3);
     alpha=compute_alpha(u, params);
 
+    /*
     // to avoid errors building up in imaginary part
     for(kx=-params.K;kx<=params.K;kx++){
       for(ky=-params.K;ky<=params.K;ky++){
 	u[KLOOKUP(kx,ky,params.S)]=__real__ u[KLOOKUP(kx,ky,params.S)];
       }
     }
+    */
 
     // running average
     if(t>0){
       avg=avg-(avg-alpha)/t;
     }
 
-    if(t>0 && t%1000==0){
+    if(t>0 && t%1==0){
       fprintf(stderr,"%8d % .8e % .8e % .8e % .8e\n",t, __real__ avg, __imag__ avg, __real__ alpha, __imag__ alpha);
       printf("%8d % .8e % .8e % .8e % .8e\n",t, __real__ avg, __imag__ avg, __real__ alpha, __imag__ alpha);
     }
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;
       }
     }
   }