Values that seem to work

1 files changed, 46 insertions, 15 deletions

diff --git a/src/main.c b/src/main.c
index e92d12f..4b68f92 100644
--- a/src/main.c
+++ b/src/main.c
@@ -72,15 +72,15 @@ int read_args(int argc, const char* argv[], ns_params* params, unsigned int* nst
 
   // defaults
   /*
-  params->h=6.103515625e-05;
-  params->K=16;
-  *nsteps=16777216;
-  params->nu=4.9632717887631524e-05;
-  */
   params->K=16;
   params->h=1e-3/(2*params->K+1);
   *nsteps=10000000;
   params->nu=1./1024/(2*params->K+1);
+  */
+  params->K=3;
+  params->h=0.0001220703125;
+  *nsteps=10000000;
+  params->nu=0.00048828125;
 
   // loop over arguments
   for(i=1;i<argc;i++){
@@ -190,6 +190,7 @@ int enstrophy(ns_params params, unsigned int Nsteps){
   unsigned int t;
   int kx,ky;
   fft_vects fft_vects;
+  double rescale;
 
   // sizes
   params.S=2*params.K+1;
@@ -203,25 +204,55 @@ int enstrophy(ns_params params, unsigned int Nsteps){
   tmp2=calloc(sizeof(_Complex double),params.S*params.S);
   tmp3=calloc(sizeof(_Complex double),params.S*params.S);
 
+  srand(17);
+
   // initial value
+  for(ky=0;ky<=params.K;ky++){
+    u[KLOOKUP(0,ky,params.S)]=(-RAND_MAX*0.5+rand())*1.0/RAND_MAX+(-RAND_MAX*0.5+rand())*1.0/RAND_MAX*I;
+  }
+  for(kx=1;kx<=params.K;kx++){
+    for(ky=-params.K;ky<=params.K;ky++){
+      u[KLOOKUP(kx,ky,params.S)]=(-RAND_MAX*0.5+rand())*1.0/RAND_MAX+(-RAND_MAX*0.5+rand())*1.0/RAND_MAX*I;
+    }
+  }
+  for(ky=-params.K;ky<=-1;ky++){
+    u[KLOOKUP(0,ky,params.S)]=conj(u[KLOOKUP(0,-ky,params.S)]);
+  }
+  for(kx=-params.K;kx<=-1;kx++){
+    for(ky=-params.K;ky<=params.K;ky++){
+      u[KLOOKUP(kx,ky,params.S)]=conj(u[KLOOKUP(-kx,-ky,params.S)]);
+    }
+  }
+  rescale=0;
   for(kx=-params.K;kx<=params.K;kx++){
     for(ky=-params.K;ky<=params.K;ky++){
-      //u[KLOOKUP(kx,ky,params.S)]=kx*kx*ky*ky*exp(-(kx*kx+ky*ky));
-      if((kx==1 && ky==0) || (kx==-1 && ky==0)){
-	u[KLOOKUP(kx,ky,params.S)]=1;
-      }
-      else{
-	u[KLOOKUP(kx,ky,params.S)]=0;
-      }
+      rescale=rescale+((__real__ u[KLOOKUP(kx,ky,params.S)])*(__real__ u[KLOOKUP(kx,ky,params.S)])+(__imag__ u[KLOOKUP(kx,ky,params.S)])*(__imag__ u[KLOOKUP(kx,ky,params.S)]))*(kx*kx+ky*ky);
+    }
+  }
+  for(kx=-params.K;kx<=params.K;kx++){
+    for(ky=-params.K;ky<=params.K;ky++){
+      u[KLOOKUP(kx,ky,params.S)]=u[KLOOKUP(kx,ky,params.S)]*sqrt(155.1/rescale);
+    }
+  }
+
+  /*
+  for(kx=-params.K;kx<=params.K;kx++){
+    for(ky=-params.K;ky<=params.K;ky++){
+      printf("%d %d % .8e % .8e\n",kx,ky, __real__ u[KLOOKUP(kx,ky,params.S)], __imag__ u[KLOOKUP(kx,ky,params.S)]);
     }
   }
+  */
+
 
   // driving force
   for(kx=-params.K;kx<=params.K;kx++){
     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.0*params.K;
+      if(kx==2 && ky==-1){
+	params.g[KLOOKUP(kx,ky,params.S)]=0.5+sqrt(3)/2*I;
+      }
+      else if(kx==-2 && ky==1){
+	params.g[KLOOKUP(kx,ky,params.S)]=0.5-sqrt(3)/2*I;
       }
       else{
 	params.g[KLOOKUP(kx,ky,params.S)]=0;
@@ -260,7 +291,7 @@ int enstrophy(ns_params params, unsigned int Nsteps){
       avg=avg-(avg-alpha)/t;
     }
 
-    if(t>0 && t%1000==0){
+    if(t>0 && t%8192==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);
     }