Only store u[kx,ky] with kx>=0

1 files changed, 58 insertions, 44 deletions

diff --git a/src/navier-stokes.c b/src/navier-stokes.c
index 0fc541e..19e2f69 100644
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -58,10 +58,10 @@ int uk(
       for(kx=-K1;kx<=K1;kx++){
 	for (ky=-K2;ky<=K2;ky++){
 	  if (kx*kx+ky*ky<=1){
-	    fprintf(stderr,"% .8e % .8e ",__real__ u[klookup(kx,ky,2*K1+1,2*K2+1)], __imag__ u[klookup(kx,ky,2*K1+1,2*K2+1)]);
+	    fprintf(stderr,"% .8e % .8e ",__real__ getval_sym(u,kx,ky,K2), __imag__ getval_sym(u, kx,ky,K2));
 	    
 	  }
-	  printf("% .15e % .15e ",__real__ u[klookup(kx,ky,2*K1+1,2*K2+1)], __imag__ u[klookup(kx,ky,2*K1+1,2*K2+1)]);
+	  printf("% .15e % .15e ",__real__ getval_sym(u, kx,ky,K2), __imag__ getval_sym(u, kx,ky,K2));
 	}
       }
       fprintf(stderr,"\n");
@@ -204,12 +204,12 @@ int  ns_init_tmps(
   unsigned int nthreads
 ){
   // velocity field
-  *u=calloc(sizeof(_Complex double),(2*K1+1)*(2*K2+1));
+  *u=calloc(sizeof(_Complex double),(K1+1)*(2*K2+1));
 
   // allocate tmp vectors for computation
-  *tmp1=calloc(sizeof(_Complex double),(2*K1+1)*(2*K2+1));
-  *tmp2=calloc(sizeof(_Complex double),(2*K1+1)*(2*K2+1));
-  *tmp3=calloc(sizeof(_Complex double),(2*K1+1)*(2*K2+1));
+  *tmp1=calloc(sizeof(_Complex double),(K1+1)*(2*K2+1));
+  *tmp2=calloc(sizeof(_Complex double),(K1+1)*(2*K2+1));
+  *tmp3=calloc(sizeof(_Complex double),(K1+1)*(2*K2+1));
 
   // init threads
   fftw_init_threads();
@@ -266,7 +266,7 @@ int copy_u(
 ){
   int i;
 
-  for(i=0;i<(2*K1+1)*(2*K2+1);i++){
+  for(i=0;i<(K1+1)*(2*K2+1);i++){
     u[i]=u0[i];
   }
 
@@ -297,61 +297,54 @@ int ns_step(
   // k1
   ns_rhs(tmp1, u, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
   // add to output
-  for(kx=-K1;kx<=K1;kx++){
+  for(kx=0;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
-      tmp3[klookup(kx,ky,2*K1+1,2*K2+1)]=u[klookup(kx,ky,2*K1+1,2*K2+1)]+delta/6*tmp1[klookup(kx,ky,2*K1+1,2*K2+1)];
+      tmp3[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta/6*tmp1[klookup_sym(kx,ky,K2)];
     }
   }
 
   // u+h*k1/2
-  for(kx=-K1;kx<=K1;kx++){
+  for(kx=0;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
-      tmp2[klookup(kx,ky,2*K1+1,2*K2+1)]=u[klookup(kx,ky,2*K1+1,2*K2+1)]+delta/2*tmp1[klookup(kx,ky,2*K1+1,2*K2+1)];
+      tmp2[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta/2*tmp1[klookup_sym(kx,ky,K2)];
     }
   }
   // k2
   ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
   // add to output
-  for(kx=-K1;kx<=K1;kx++){
+  for(kx=0;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
-      tmp3[klookup(kx,ky,2*K1+1,2*K2+1)]+=delta/3*tmp1[klookup(kx,ky,2*K1+1,2*K2+1)];
+      tmp3[klookup_sym(kx,ky,K2)]+=delta/3*tmp1[klookup_sym(kx,ky,K2)];
     }
   }
 
   // u+h*k2/2
-  for(kx=-K1;kx<=K1;kx++){
+  for(kx=0;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
-      tmp2[klookup(kx,ky,2*K1+1,2*K2+1)]=u[klookup(kx,ky,2*K1+1,2*K2+1)]+delta/2*tmp1[klookup(kx,ky,2*K1+1,2*K2+1)];
+      tmp2[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta/2*tmp1[klookup_sym(kx,ky,K2)];
     }
   }
   // k3
   ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
   // add to output
-  for(kx=-K1;kx<=K1;kx++){
+  for(kx=0;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
-      tmp3[klookup(kx,ky,2*K1+1,2*K2+1)]+=delta/3*tmp1[klookup(kx,ky,2*K1+1,2*K2+1)];
+      tmp3[klookup_sym(kx,ky,K2)]+=delta/3*tmp1[klookup_sym(kx,ky,K2)];
     }
   }
 
   // u+h*k3
-  for(kx=-K1;kx<=K1;kx++){
+  for(kx=0;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
-      tmp2[klookup(kx,ky,2*K1+1,2*K2+1)]=u[klookup(kx,ky,2*K1+1,2*K2+1)]+delta*tmp1[klookup(kx,ky,2*K1+1,2*K2+1)];
+      tmp2[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta*tmp1[klookup_sym(kx,ky,K2)];
     }
   }
   // k4
   ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
   // add to output
-  for(kx=-K1;kx<=K1;kx++){
+  for(kx=0;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
-      u[klookup(kx,ky,2*K1+1,2*K2+1)]=tmp3[klookup(kx,ky,2*K1+1,2*K2+1)]+delta/6*tmp1[klookup(kx,ky,2*K1+1,2*K2+1)];
-    }
-  }
-
-  // enforce u(-k)=u(k)^*
-  for(kx=-K1;kx<=K1;kx++){
-    for(ky=-K2;ky<=K2;ky++){
-      u[klookup(kx,ky,2*K1+1,2*K2+1)]=(u[klookup(kx,ky,2*K1+1,2*K2+1)]+conj(u[klookup(-kx,-ky,2*K1+1,2*K2+1)]))/2;
+      u[klookup_sym(kx,ky,K2)]=tmp3[klookup_sym(kx,ky,K2)]+delta/6*tmp1[klookup_sym(kx,ky,K2)];
     }
   }
 
@@ -397,13 +390,13 @@ int ns_rhs(
   printf("% .15e\n",sqrt(cmp));
   */
 
-  for(i=0; i<(2*K1+1)*(2*K2+1); i++){
+  for(i=0; i<(K1+1)*(2*K2+1); i++){
     out[i]=0;
   }
-  for(kx=-K1;kx<=K1;kx++){
+  for(kx=0;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
       if(kx!=0 || ky!=0){
-        out[klookup(kx,ky,2*K1+1,2*K2+1)]=-4*M_PI*M_PI/L/L*alpha*(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]+g[klookup(kx,ky,2*K1+1,2*K2+1)]+4*M_PI*M_PI/L/L/sqrt(kx*kx+ky*ky)*ifft.fft[klookup(kx,ky,N1,N2)];
+        out[klookup_sym(kx,ky,K2)]=-4*M_PI*M_PI/L/L*alpha*(kx*kx+ky*ky)*u[klookup_sym(kx,ky,K2)]+g[klookup_sym(kx,ky,K2)]+4*M_PI*M_PI/L/L/sqrt(kx*kx+ky*ky)*ifft.fft[klookup(kx,ky,N1,N2)];
       }
     }
   }
@@ -436,8 +429,8 @@ int ns_T(
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
       if(kx!=0 || ky!=0){
-        fft1.fft[klookup(kx,ky,N1,N2)]=kx/sqrt(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]/N1;
-        fft2.fft[klookup(kx,ky,N1,N2)]=ky*sqrt(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]/N2;
+        fft1.fft[klookup(kx,ky,N1,N2)]=kx/sqrt(kx*kx+ky*ky)*getval_sym(u, kx,ky,K2)/N1;
+        fft2.fft[klookup(kx,ky,N1,N2)]=ky*sqrt(kx*kx+ky*ky)*getval_sym(u, kx,ky,K2)/N2;
       }
     }
   }
@@ -460,8 +453,8 @@ int ns_T(
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
       if(kx!=0 || ky!=0){
-        fft1.fft[klookup(kx,ky,N1,N2)]=ky/sqrt(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]/N1;
-        fft2.fft[klookup(kx,ky,N1,N2)]=kx*sqrt(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]/N2;
+        fft1.fft[klookup(kx,ky,N1,N2)]=ky/sqrt(kx*kx+ky*ky)*getval_sym(u, kx,ky,K2)/N1;
+        fft2.fft[klookup(kx,ky,N1,N2)]=kx*sqrt(kx*kx+ky*ky)*getval_sym(u, kx,ky,K2)/N2;
       }
     }
   }
@@ -477,14 +470,12 @@ int ns_T(
   // inverse fft
   fftw_execute(ifft.fft_plan);
 
-  /*
   // enforce T(u,-k)=T(u,k)^*
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
       ifft.fft[klookup(kx,ky,N1,N2)]=(ifft.fft[klookup(kx,ky,N1,N2)]+conj(ifft.fft[klookup(-kx,-ky,N1,N2)]))/2;
     }
   }
-  */
 
   return(0);
 }
@@ -519,7 +510,7 @@ int ns_T_nofft(
 
 	  // cutoff in q
 	  if(qx>=-K1 && qx<=K1 && qy>=-K2 && qy<=K2 && qx*qx+qy*qy>0 && px*px+py*py>0){
-	    out[klookup(kx,ky,N1,N2)]+=(-qx*py+qy*px)*sqrt(qx*qx+qy*qy)/sqrt(px*px+py*py)*u[klookup(px,py,2*K1+1,2*K2+1)]*u[klookup(qx,qy,2*K1+1,2*K2+1)];
+	    out[klookup(kx,ky,N1,N2)]+=(-qx*py+qy*px)*sqrt(qx*qx+qy*qy)/sqrt(px*px+py*py)*getval_sym(u, px,py,K2)*getval_sym(u, qx,qy,K2);
 	  }
 	}
       }
@@ -556,8 +547,8 @@ double compute_alpha(
 
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
-      num+=(L*L/4/M_PI/M_PI*(kx*kx+ky*ky)*g[klookup(kx,ky,2*K1+1,2*K2+1)]+sqrt(kx*kx+ky*ky)*T[klookup(kx,ky,N1,N2)])*conj(u[klookup(kx,ky,2*K1+1,2*K2+1)]);
-      denom+=__real__ (kx*kx+ky*ky)*(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]*conj(u[klookup(kx,ky,2*K1+1,2*K2+1)]);
+      num+=(L*L/4/M_PI/M_PI*(kx*kx+ky*ky)*getval_sym(g, kx,ky,K2)+sqrt(kx*kx+ky*ky)*T[klookup(kx,ky,N1,N2)])*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));
     }
   }
 
@@ -579,8 +570,8 @@ double compute_alpha(
 
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
-      denom+=(kx*kx+ky*ky)*(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]*conj(u[klookup(kx,ky,2*K1+1,2*K2+1)])*(1+(ky!=0?kx*kx/ky/ky:0));
-      num+=(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]*conj(g[klookup(kx,ky,2*K1+1,2*K2+1)])*(1+(ky!=0?kx*kx/ky/ky:0));
+      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));
     }
   }
 
@@ -597,7 +588,7 @@ double compute_energy(
   double out=0.;
   for(kx=-K1;kx<=K1;kx++){
     for (ky=-K2;ky<=K2;ky++){
-      out+=__real__ (u[klookup(kx,ky,2*K1+1,2*K2+1)]*conj(u[klookup(kx,ky,2*K1+1,2*K2+1)]));
+      out+=__real__ (getval_sym(u, kx,ky,K2)*conj(getval_sym(u, kx,ky,K2)));
     }
   }
   return out;
@@ -614,7 +605,7 @@ double compute_enstrophy(
   double out=0.;
   for(kx=-K1;kx<=K1;kx++){
     for (ky=-K2;ky<=K2;ky++){
-      out+=__real__ (4*M_PI*M_PI/L/L*(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]*conj(u[klookup(kx,ky,2*K1+1,2*K2+1)]));
+      out+=__real__ (4*M_PI*M_PI/L/L*(kx*kx+ky*ky)*getval_sym(u, kx,ky,K2)*conj(getval_sym(u, kx,ky,K2)));
     }
   }
   return out;
@@ -630,3 +621,26 @@ int klookup(
 ){
   return (kx>=0 ? kx : S1+kx)*S2 + (ky>=0 ? ky : S2+ky);
 }
+
+// get index for kx,ky in array of size K1,K2 in which only the terms with kx>=0 are stored
+int klookup_sym(
+  int kx,
+  int ky,
+  int K2
+){
+  if (kx<0) {
+    fprintf(stderr, "bug!: attempting to access a symmetrized vector at kx<0\n Contact Ian at ian.jauslin@rutgers.edu!\n");
+    exit(-1);
+  }
+  return kx*(2*K2+1) + (ky>=0 ? ky : (2*K2+1)+ky);
+}
+
+// get u_{kx,ky} from a vector u in which only the values for kx>=0 are stored, assuming u_{-k}=u_k^*
+_Complex double getval_sym(
+  _Complex double* u,
+  int kx,
+  int ky,
+  int K2
+){
+  return (kx>=0 ? u[klookup_sym(kx,ky,K2)] : conj(u[klookup_sym(-kx,-ky,K2)]));
+}