diff options
| -rw-r--r-- | src/main.c | 61 |
1 files changed, 46 insertions, 15 deletions
@@ -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); } |