Remove terms with kx=0 and ky<=0

5 files changed, 50 insertions, 37 deletions

diff --git a/src/driving.c b/src/driving.c
index 7920d45..527a37d 100644
--- a/src/driving.c
+++ b/src/driving.c
@@ -11,7 +11,7 @@ int g_test(
 ){
   int kx,ky;
   for(kx=0;kx<=K1;kx++){
-    for (ky=-K2;ky<=K2;ky++){
+    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
       if(kx==2 && ky==-1){
 	g[klookup_sym(kx,ky,K2)]=0.5+sqrt(3)/2*I;
       }
@@ -30,7 +30,7 @@ int g_zero(
   int K2
 ){
   int i;
-  for(i=0;i<(K1+1)*(2*K2+1);i++){
+  for(i=0;i<K1*(2*K2+1)+K2;i++){
     g[i]=0.;
   }
 
diff --git a/src/init.c b/src/init.c
index cce97d2..cc482da 100644
--- a/src/init.c
+++ b/src/init.c
@@ -22,12 +22,10 @@ int init_random (
 
   // random init (set half, then the other half are the conjugates)
   for(kx=0;kx<=K1;kx++){
-    for(ky=-K2;ky<=K2;ky++){
-      if (kx!=0 || ky>0){
-	x=-0.5+((double) rand())/RAND_MAX;
-	y=-0.5+((double) rand())/RAND_MAX;
-	u0[klookup_sym(kx,ky,K2)]=x+y*I;
-      }
+    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+      x=-0.5+((double) rand())/RAND_MAX;
+      y=-0.5+((double) rand())/RAND_MAX;
+      u0[klookup_sym(kx,ky,K2)]=x+y*I;
     }
   }
 
@@ -39,7 +37,7 @@ int init_random (
     rescale=compute_enstrophy(u0, K1, K2, L);
   }
   for(kx=0;kx<=K1;kx++){
-    for(ky=-K2;ky<=K2;ky++){
+    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
       u0[klookup_sym(kx,ky,K2)]=u0[klookup_sym(kx,ky,K2)]*sqrt(init_en/rescale);
     }
   }
@@ -60,7 +58,7 @@ int init_gaussian (
   double rescale;
 
   for(kx=0;kx<=K1;kx++){
-    for(ky=-K2;ky<=K2;ky++){
+    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
       u0[klookup_sym(kx,ky,K2)]=(kx*kx+ky*ky)*exp(-(kx*kx+ky*ky));
     }
   }
@@ -73,7 +71,7 @@ int init_gaussian (
     rescale=compute_enstrophy(u0, K1, K2, L);
   }
   for(kx=0;kx<=K1;kx++){
-    for(ky=-K2;ky<=K2;ky++){
+    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
       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 e8afbb1..c01f0c9 100644
--- a/src/io.c
+++ b/src/io.c
@@ -14,7 +14,7 @@ int write_vec(_Complex double* vec, int K1, int K2, FILE* file){
   }
 
   for(kx=0;kx<=K1;kx++){
-    for (ky=-K2;ky<=K2;ky++){
+    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
       fprintf(file,"% 3d % 3d % .15e % .15e\n",kx,ky,__real__ vec[klookup_sym(kx,ky,K2)],__imag__ vec[klookup_sym(kx,ky,K2)]);
     }
   }
@@ -29,7 +29,7 @@ int write_vec_bin(_Complex double* vec, int K1, int K2, FILE* file){
     return 0;
   }
 
-  fwrite(vec, sizeof(_Complex double), (K1+1)*(2*K2+1), file);
+  fwrite(vec, sizeof(_Complex double), K1*(2*K2+1)+K2, file);
 
   return 0;
 }
@@ -83,7 +83,10 @@ int read_vec(_Complex double* out, int K1, int K2, FILE* file){
 	    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);
+	    fprintf(stderr, "warning: reading line %d: kx should be >=0, got %d, skipping\n", counter, kx);
+	  }
+	  else if (kx==0 && ky<=0){
+	    fprintf(stderr, "warning: reading line %d: if kx==0 then ky should be >0, got kx=%d ky=%d, skipping\n", counter, kx, ky);
 	  }
 	  else{
 	    // set output
@@ -164,7 +167,7 @@ int read_vec_bin(_Complex double* out, int K1, int K2, FILE* file){
     }
   }
 
-  fread(out, sizeof(_Complex double), (K1+1)*(2*K2+1), file);
+  fread(out, sizeof(_Complex double), K1*(2*K2+1)+K2, file);
 
   return 0;
 }
diff --git a/src/main.c b/src/main.c
index 7b3cdf2..d7d30bb 100644
--- a/src/main.c
+++ b/src/main.c
@@ -635,7 +635,7 @@ int set_parameter(
 _Complex double* set_driving(
   nstrophy_parameters parameters
 ){
-  _Complex double* g=calloc(sizeof(_Complex double),(2*parameters.K1+1)*(2*parameters.K2+1));
+  _Complex double* g=calloc(sizeof(_Complex double),parameters.K1*(2*parameters.K2+1)+parameters.K2);
 
   switch(parameters.driving){
   case DRIVING_ZERO:
@@ -664,7 +664,7 @@ _Complex double* set_driving(
 _Complex double* set_init(
   nstrophy_parameters parameters
 ){
-  _Complex double* u0=calloc(sizeof(_Complex double),(2*parameters.K1+1)*(2*parameters.K2+1));
+  _Complex double* u0=calloc(sizeof(_Complex double),parameters.K1*(2*parameters.K2+1)+parameters.K2);
 
   switch(parameters.init){
   case INIT_RANDOM:
diff --git a/src/navier-stokes.c b/src/navier-stokes.c
index 153631b..51d198c 100644
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -241,12 +241,12 @@ int  ns_init_tmps(
   unsigned int nthreads
 ){
   // velocity field
-  *u=calloc(sizeof(_Complex double),(K1+1)*(2*K2+1));
+  *u=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
 
   // allocate tmp vectors for computation
-  *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));
+  *tmp1=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
+  *tmp2=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
+  *tmp3=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
 
   // init threads
   fftw_init_threads();
@@ -303,7 +303,7 @@ int copy_u(
 ){
   int i;
 
-  for(i=0;i<(K1+1)*(2*K2+1);i++){
+  for(i=0;i<K1*(2*K2+1)+K2;i++){
     u[i]=u0[i];
   }
 
@@ -335,14 +335,14 @@ int ns_step(
   ns_rhs(tmp1, u, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
   // add to output
   for(kx=0;kx<=K1;kx++){
-    for(ky=-K2;ky<=K2;ky++){
+    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
       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=0;kx<=K1;kx++){
-    for(ky=-K2;ky<=K2;ky++){
+    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
       tmp2[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta/2*tmp1[klookup_sym(kx,ky,K2)];
     }
   }
@@ -350,14 +350,14 @@ int ns_step(
   ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
   // add to output
   for(kx=0;kx<=K1;kx++){
-    for(ky=-K2;ky<=K2;ky++){
+    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
       tmp3[klookup_sym(kx,ky,K2)]+=delta/3*tmp1[klookup_sym(kx,ky,K2)];
     }
   }
 
   // u+h*k2/2
   for(kx=0;kx<=K1;kx++){
-    for(ky=-K2;ky<=K2;ky++){
+    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
       tmp2[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta/2*tmp1[klookup_sym(kx,ky,K2)];
     }
   }
@@ -365,14 +365,14 @@ int ns_step(
   ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
   // add to output
   for(kx=0;kx<=K1;kx++){
-    for(ky=-K2;ky<=K2;ky++){
+    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
       tmp3[klookup_sym(kx,ky,K2)]+=delta/3*tmp1[klookup_sym(kx,ky,K2)];
     }
   }
 
   // u+h*k3
   for(kx=0;kx<=K1;kx++){
-    for(ky=-K2;ky<=K2;ky++){
+    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
       tmp2[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta*tmp1[klookup_sym(kx,ky,K2)];
     }
   }
@@ -380,7 +380,7 @@ int ns_step(
   ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
   // add to output
   for(kx=0;kx<=K1;kx++){
-    for(ky=-K2;ky<=K2;ky++){
+    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
       u[klookup_sym(kx,ky,K2)]=tmp3[klookup_sym(kx,ky,K2)]+delta/6*tmp1[klookup_sym(kx,ky,K2)];
     }
   }
@@ -417,14 +417,12 @@ int ns_rhs(
     alpha=compute_alpha(u,K1,K2,g,L);
   }
 
-  for(i=0; i<(K1+1)*(2*K2+1); i++){
+  for(i=0; i<K1*(2*K2+1)+K2; i++){
     out[i]=0;
   }
   for(kx=0;kx<=K1;kx++){
-    for(ky=-K2;ky<=K2;ky++){
-      if(kx!=0 || ky!=0){
-        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)];
-      }
+    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+      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)];
     }
   }
 
@@ -610,7 +608,7 @@ 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
+// get index for kx,ky in array of size K1,K2 in which only the terms with kx>=0 and if kx=0, ky>0 are stored
 int klookup_sym(
   int kx,
   int ky,
@@ -620,7 +618,14 @@ int klookup_sym(
     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);
+  if (kx==0) {
+    if (ky<=0){
+      fprintf(stderr, "bug!: attempting to access a symmetrized vector at kx=0 and ky<=0\n Contact Ian at ian.jauslin@rutgers.edu!\n");
+      exit(-1);
+    }
+    return ky-1;
+  }
+  return K2+(kx-1)*(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^*
@@ -630,5 +635,12 @@ _Complex double getval_sym(
   int ky,
   int K2
 ){
-  return (kx>=0 ? u[klookup_sym(kx,ky,K2)] : conj(u[klookup_sym(-kx,-ky,K2)]));
+  if(kx>0 || (kx==0 && ky>0)){
+    return u[klookup_sym(kx,ky,K2)];
+  }
+  else if(kx==0 && ky==0){
+    return 0;
+  } else {
+    return conj(u[klookup_sym(-kx,-ky,K2)]);
+  }
 }