Only use T term in alpha for ns_step computation

1 files changed, 14 insertions, 23 deletions

diff --git a/src/navier-stokes.c b/src/navier-stokes.c
index ab56a7e..78207e0 100644
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -128,11 +128,8 @@ int eea(
   for(t=starting_time;nsteps==0 || t<starting_time+nsteps;t++){
     ns_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
 
-    // convolution term
-    ns_T(u,K1,K2,N1,N2,fft1,fft2,ifft);
-
     energy=compute_energy(u, K1, K2);
-    alpha=compute_alpha(u, K1, K2, N1, N2, g, L, ifft.fft);
+    alpha=compute_alpha_fast(u, K1, K2, g, L);
     enstrophy=compute_enstrophy(u, K1, K2, L);
     
     // running average
@@ -434,16 +431,6 @@ int ns_rhs(
     alpha=compute_alpha(u,K1,K2,N1,N2,g,L,ifft.fft);
   }
 
-  /*
-  // compare convolution term (store result in fft1.fft)
-  ns_T_nofft(fft1.fft, u, K1, K2, N1, N2);
-  double cmp=0.;
-  for(i=0;i<N1*N2;i++){
-    cmp+=(ifft.fft[i]-fft1.fft[i])*(ifft.fft[i]-fft1.fft[i]);
-  }
-  printf("% .15e\n",sqrt(cmp));
-  */
-
   for(i=0; i<(K1+1)*(2*K2+1); i++){
     out[i]=0;
   }
@@ -609,28 +596,32 @@ double compute_alpha(
   return __real__ num/denom;
 }
 
-
-/*
 // compute alpha
-double compute_alpha(
+// do not include the term involving T, which, in the limit K->infty, vanishes
+double compute_alpha_fast(
   _Complex double* u,
   int K1,
   int K2,
-  _Complex double* g
+  _Complex double* g,
+  double L
 ){
   _Complex double num=0;
-  _Complex double denom=0;
+  double denom=0;
   int kx,ky;
 
+  num=0.;
+  denom=0.;
+
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
-      denom+=(kx*kx+ky*ky)*(kx*kx+ky*ky)*getval_sym(u, kx,ky,K2)*conj(getval_sym(u, kx,ky,K2))*(1+(ky!=0?kx*kx/ky/ky:0));
-      num+=(kx*kx+ky*ky)*getval_sym(u, kx,ky,K2)*conj(g[getval(kx,ky,K1,K2)])*(1+(ky!=0?kx*kx/ky/ky:0));
+      num+=L*L/4/M_PI/M_PI*(kx*kx+ky*ky)*getval_sym(g, kx,ky,K2)*conj(getval_sym(u, kx,ky,K2));
+      denom+=__real__ (kx*kx+ky*ky)*(kx*kx+ky*ky)*getval_sym(u, kx,ky,K2)*conj(getval_sym(u, kx,ky,K2));
     }
   }
 
-  return __real__ (num/denom);
-}*/
+  return __real__ num/denom;
+}
+
 
 // compute energy
 double compute_energy(