Remove old ways of computing convolution

1 files changed, 0 insertions, 136 deletions

diff --git a/src/navier-stokes.c b/src/navier-stokes.c
index 3def963..3dd7484 100644
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -3,8 +3,6 @@
 
 #define M_PI 3.14159265358979323846
 
-#define CHK 1
-
 // next time step for Irreversible Navier-Stokes equation
 int ins_step(_Complex double* u, ns_params params, fft_vects vects, _Complex double* tmp1, _Complex double* tmp2, _Complex double* tmp3){
   int kx,ky;
@@ -70,7 +68,6 @@ 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;
 
-#if CHK==0
   // F(u/|p|)*F(q1*q2*u/|q|)
   // init to 0
   for(kx=0; kx<params.N*params.N; kx++){
@@ -144,140 +141,7 @@ int ins_rhs(_Complex double* out, _Complex double* u, ns_params params, fft_vect
       }
     }
   }
-#elif CHK == 1
-  // F(-p2/|p|*u)*F(q1*|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.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)]=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.fft1[KLOOKUP(ky,kx,params.N)]*vects.fft2[KLOOKUP(ky,kx,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/params.S*(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)]+params.g[KLOOKUP(kx,ky,params.S)]+2*M_PI/sqrt(kx*kx+ky*ky)/params.S*vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N;
-      }
-    }
-  }
-
-#elif CHK==2
-  // F(u)*F(q1*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++){
-      // u_k
-      vects.fft1[KLOOKUP(kx,ky,params.N)]=u[KLOOKUP(kx,ky,params.S)];
-      // 2i\pi k_x u_k
-      vects.fft2[KLOOKUP(kx,ky,params.N)]=2*M_PI*I*kx*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)];
-    }
-  }
-
-  // F(p1/p2*u)*F(q2*u)
-  // init to 0
-  for(kx=0; kx<params.N*params.N; kx++){
-    vects.fft1[kx]=0;
-    vects.fft2[kx]=0;
-  }
-  // fill modes
-  for(kx=-params.K;kx<=params.K;kx++){
-    for(ky=-params.K;ky<=params.K;ky++){
-      // k_x/k_y u_k
-      if(ky!=0){
-	vects.fft1[KLOOKUP(kx,ky,params.N)]=kx/ky*u[KLOOKUP(kx,ky,params.S)];
-      }
-      // 2i\pi k_y u_k
-      vects.fft2[KLOOKUP(kx,ky,params.N)]=2*M_PI*I*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)];
-    }
-  }
-
-  // inverse fft
-  fftw_execute(vects.invfft_plan);
-
-  /*
-  // check: convolution instead of fft
-  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++){
-      for(px=-params.K;px<=params.K;px++){
-	for(py=-params.K;py<=params.K;py++){
-	  if(kx-px<=params.K && kx-px>=-params.K && ky-py<=params.K && ky-py>=-params.K){
-	    out[KLOOKUP(kx,ky,params.S)]+=2*M_PI*I*(u[KLOOKUP(px,py,params.S)]*(kx-px)*u[KLOOKUP(kx-px,ky-py,params.S)]-(py==0?0:px/py*u[KLOOKUP(px,py,params.S)]*(ky-py)*u[KLOOKUP(kx-px,ky-py,params.S)]));
-	  }
-	}
-      }
-      dd=(__real__ vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N-__real__ out[KLOOKUP(kx,ky,params.S)])*(__real__ vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N-__real__ out[KLOOKUP(kx,ky,params.S)])+(__imag__ vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N-__imag__ out[KLOOKUP(kx,ky,params.S)])*(__imag__ vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N-__imag__ out[KLOOKUP(kx,ky,params.S)]);
-      if(dd>1e-25){
-	printf("%d %d % .8e\n",kx,ky, dd);
-      }
-    }
-  }
-  */
-
-
-  // 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++){
-      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)]+vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N;
-    }
-  }
 
-#endif
   return(0);
 }