Simpler expression for fft term

docs/nstrophy_doc/nstrophy_doc.tex

@@ -49,16 +49,23 @@ in terms of which, multiplying both sides of the equation by $\frac{k^\perp}{|k|
 \begin{equation}
   \partial_t\hat \varphi_k=
   -4\pi^2\nu k^2\hat \varphi_k+\frac{k^\perp}{2i\pi|k|}\cdot\hat g_k
-  -\frac{2i\pi}{|k|}\sum_{\displaystyle\mathop{\scriptstyle p,q\in\mathbb Z^2}_{p+q=k}}
+  +\frac{4\pi^2}{|k|}\sum_{\displaystyle\mathop{\scriptstyle p,q\in\mathbb Z^2}_{p+q=k}}
   (q\cdot p^\perp)(k^\perp\cdot q^\perp)\hat\varphi_p\hat\varphi_q
+  .
   \label{ins_k}
 \end{equation}
-which, since $q\cdot p^\perp=q\cdot(k^\perp-q^\perp)=q\cdot k^\perp$, is
+Furthermore
+\begin{equation}
+  (q\cdot p^\perp)(k^\perp\cdot q^\perp)
+  =
+  (q\cdot p^\perp)(q^2+p\cdot q)
+\end{equation}
+and $q\cdot p^\perp$ is antisymmetric under exchange of $q$ and $p$. Therefore,
 \begin{equation}
   \partial_t\hat \varphi_k=
   -4\pi^2\nu k^2\hat \varphi_k+\frac{k^\perp}{2i\pi|k|}\cdot\hat g_k
   +\frac{4\pi^2}{|k|}\sum_{\displaystyle\mathop{\scriptstyle p,q\in\mathbb Z^2}_{p+q=k}}
-  (q\cdot k^\perp)(k\cdot q)\frac{\hat\varphi_p\hat\varphi_q}{|p||q|}
+  \frac{(q\cdot p^\perp)|q|}{|p|}\hat\varphi_p\hat\varphi_q
   .
   \label{ins_k}
 \end{equation}
@@ -73,8 +80,7 @@ We truncate the Fourier modes and assume that $\hat\varphi_k=0$ if $|k_1|>K_1$ o
 \begin{equation}
   T(\hat\varphi,k):=
   \sum_{\displaystyle\mathop{\scriptstyle p,q\in\mathbb Z^2}_{p+q=k}}
-  \frac{\hat\varphi_p}{|p|}
+  \frac{(q\cdot p^\perp)|q|}{|p|}\hat\varphi_q\hat\varphi_p
-  (q\cdot k^\perp)(k\cdot q)\frac{\hat\varphi_q}{|q|}
 \end{equation}
 using a fast Fourier transform, defined as
 \begin{equation}
@@ -98,8 +104,7 @@ in which $\mathcal F^*$ is defined like $\mathcal F$ but with the opposite phase
   \sum_{p,q\in\mathcal K}
   \frac1{N_1N_2}
   \sum_{n\in\mathcal N}e^{-\frac{2i\pi}{N_1}n_1(p_1+q_1-k_1)-\frac{2i\pi}{N_2}n_2(p_2+q_2-k_2)}
-  \frac{\hat\varphi_p}{|p|}
+  (q\cdot p^\perp)\frac{|q|}{|p|}\hat\varphi_q\hat\varphi_p
-  (q\cdot k^\perp)(k\cdot q)\frac{\hat\varphi_q}{|q|}
 \end{equation}
 provided
 \begin{equation}
@@ -109,29 +114,17 @@ Indeed, $\sum_{n_i=0}^{N_i}e^{-\frac{2i\pi}{N_i}n_im_i}$ vanishes unless $m_i=0\
 \begin{equation}
   T(\hat\varphi,k)
   =
+  \textstyle
   \frac1{N_1N_2}
   \mathcal F^*\left(
-    \mathcal F(|p|^{-1}\hat\varphi_p)(n)
+    \mathcal F\left(\frac{p_x\hat\varphi_p}{|p|}\right)(n)
-    \mathcal F((|q|^{-1}(q\cdot k^\perp)(k\cdot q)\hat\varphi_q)_q)(n)
+    \mathcal F\left(q_y|q|\hat\varphi_q\right)(n)
-  \right)(k)
+    -
-\end{equation}
+    \mathcal F\left(\frac{p_y\hat\varphi_p}{|p|}\right)(n)
-which we expand
+    \mathcal F\left(q_x|q|\hat\varphi_q\right)(n)
-\begin{equation}
-  T(\hat\varphi,k)
-  =
-  \textstyle
-  \frac{k_x^2-k_y^2}{N_1N_2}
-  \mathcal F^*\left(
-    \mathcal F\left(\frac{\hat\varphi_p}{|p|}\right)(n)
-    \mathcal F\left(\frac{q_xq_y}{|q|}\hat\varphi_q\right)(n)
-  \right)(k)
-  -
-  \frac{k_xk_y}{N_1N_2}
-  \mathcal F^*\left(
-    \mathcal F\left(\frac{\hat\varphi_p}{|p|}\right)(n)
-    \mathcal F\left(\frac{q_x^2-q_y^2}{|q|}\hat\varphi_q\right)(n)
   \right)(k)
 \end{equation}
+\bigskip
 
 
 \vfill

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);
     }

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.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*ky*u[KLOOKUP(kx,ky,params.S)]/sqrt(kx*kx+ky*ky);
+        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,6 +94,33 @@ 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)];
     }
   }
+
+  // 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;
+  }
+  // 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)]=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.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);
   
@@ -103,41 +131,7 @@ 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){
-        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;
+        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;
-      }
-    }
-  }
-
-  // F(u/|p|)*F((q1*q1-q2*q2)*u/|q|)
-  // init to 0
-  for(kx=0; kx<params.N*params.N; kx++){
-    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);
-      }
-    }
-  }
-  // fft
-  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);
-  
-  // write out
-  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;
       }
     }
   }