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

5 files changed, 81 insertions, 64 deletions

diff --git a/src/driving.c b/src/driving.c
index fa26b4b..c329848 100644
--- a/src/driving.c
+++ b/src/driving.c
@@ -9,16 +9,13 @@ int g_test(
   int K2
 ){
   int kx,ky;
-  for(kx=-K1;kx<=K1;kx++){
+  for(kx=0;kx<=K1;kx++){
     for (ky=-K2;ky<=K2;ky++){
       if(kx==2 && ky==-1){
-	g[klookup(kx,ky,2*K1+1,2*K2+1)]=0.5+sqrt(3)/2*I;
-      }
-      else if(kx==-2 && ky==1){
-        g[klookup(kx,ky,2*K1+1,2*K2+1)]=0.5-sqrt(3)/2*I;
+	g[klookup_sym(kx,ky,K2)]=0.5+sqrt(3)/2*I;
       }
       else{
-        g[klookup(kx,ky,2*K1+1,2*K2+1)]=0.;
+        g[klookup_sym(kx,ky,K2)]=0.;
       }
     }
   }
@@ -32,7 +29,7 @@ int g_zero(
   int K2
 ){
   int i;
-  for(i=0;i<(2*K1+1)*(2*K2+1);i++){
+  for(i=0;i<(K1+1)*(2*K2+1);i++){
     g[i]=0.;
   }
 
diff --git a/src/init.c b/src/init.c
index 4865ed9..ca693f9 100644
--- a/src/init.c
+++ b/src/init.c
@@ -26,8 +26,7 @@ int init_random (
       if (kx!=0 || ky>0){
 	x=-0.5+((double) rand())/RAND_MAX;
 	y=-0.5+((double) rand())/RAND_MAX;
-	u0[klookup(kx,ky,2*K1+1,2*K2+1)]=x+y*I;
-	u0[klookup(-kx,-ky,2*K1+1,2*K2+1)]=conj(u0[klookup(kx,ky,2*K1+1,2*K2+1)]);
+	u0[klookup_sym(kx,ky,K2)]=x+y*I;
       }
     }
   }
@@ -39,9 +38,9 @@ int init_random (
     // fix the enstrophy
     rescale=compute_enstrophy(u0, K1, K2, L);
   }
-  for(kx=-K1;kx<=K1;kx++){
+  for(kx=0;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
-      u0[klookup(kx,ky,2*K1+1,2*K2+1)]=u0[klookup(kx,ky,2*K1+1,2*K2+1)]*sqrt(init_en/rescale);
+      u0[klookup_sym(kx,ky,K2)]=u0[klookup_sym(kx,ky,K2)]*sqrt(init_en/rescale);
     }
   }
 
@@ -60,9 +59,9 @@ int init_gaussian (
   int kx,ky;
   double rescale;
 
-  for(kx=-K1;kx<=K1;kx++){
+  for(kx=0;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
-      u0[klookup(kx,ky,2*K1+1,2*K2+1)]=(kx*kx+ky*ky)*exp(-(kx*kx+ky*ky));
+      u0[klookup_sym(kx,ky,K2)]=(kx*kx+ky*ky)*exp(-(kx*kx+ky*ky));
     }
   }
 
@@ -73,9 +72,9 @@ int init_gaussian (
     // fix the enstrophy
     rescale=compute_enstrophy(u0, K1, K2, L);
   }
-  for(kx=-K1;kx<=K1;kx++){
+  for(kx=0;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
-      u0[klookup(kx,ky,2*K1+1,2*K2+1)]=u0[klookup(kx,ky,2*K1+1,2*K2+1)]*sqrt(init_en/rescale);
+      u0[klookup_sym(kx,ky,K2)]=u0[klookup_sym(kx,ky,K2)]*sqrt(init_en/rescale);
     }
   }
 
diff --git a/src/io.c b/src/io.c
index 5cdfc54..049f310 100644
--- a/src/io.c
+++ b/src/io.c
@@ -12,9 +12,9 @@ int write_u(_Complex double* u, int K1, int K2, FILE* file){
     return 0;
   }
 
-  for(kx=-K1;kx<=K1;kx++){
+  for(kx=0;kx<=K1;kx++){
     for (ky=-K2;ky<=K2;ky++){
-      fprintf(file,"%d %d % .15e % .15e\n",kx,ky,__real__ u[klookup(kx,ky,2*K1+1,2*K2+1)],__imag__ u[klookup(kx,ky,2*K1+1,2*K2+1)]);
+      fprintf(file,"%d %d % .15e % .15e\n",kx,ky,__real__ u[klookup_sym(kx,ky,K2)],__imag__ u[klookup_sym(kx,ky,K2)]);
     }
   }
 
@@ -63,15 +63,18 @@ int read_u(_Complex double* u, int K1, int K2, FILE* file){
 
 	// errors
 	if(ret!=4){
-	  fprintf(stderr, "warning: line %d does not match the input format: '%s'\n", counter, line);
+	  fprintf(stderr, "warning: line %d does not match the input format: '%s', skipping\n", counter, line);
 	}
 	else{
 	  if(kx>K1 || kx<-K1 || ky>K2 || ky<-K2){
-	    fprintf(stderr, "warning: reading line %d: kx or ky out of bounds: %d,%d\n", counter, kx, ky);
+	    fprintf(stderr, "warning: reading line %d: kx or ky out of bounds: %d,%d, skipping\n", counter, kx, ky);
+	  }
+	  else if (kx<0){
+	    fprintf(stderr, "warning: reading line %d: kx should be >=0, skipping\n", counter);
 	  }
 	  else{
 	    // set u
-	    u[klookup(kx, ky, 2*K1+1, 2*K2+1)]=r+i*I;
+	    u[klookup_sym(kx, ky, K2)]=r+i*I;
 	  }
 	}
       }
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)]));
+}
diff --git a/src/navier-stokes.h b/src/navier-stokes.h
index 582212a..d100f44 100644
--- a/src/navier-stokes.h
+++ b/src/navier-stokes.h
@@ -52,6 +52,10 @@ double compute_enstrophy( _Complex double* u, int K1, int K2, double L);
 
 // get index for kx,ky in array of size S
 int klookup( int kx, int ky, int S1, int S2);
+// 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);
+// 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);
 
 #endif