From d8380602656b54649c46bac93876b814b4d8eba7 Mon Sep 17 00:00:00 2001 From: Ian Jauslin Date: Tue, 25 Apr 2023 17:51:14 -0400 Subject: Remove terms with kx=0 and ky<=0 --- src/driving.c | 4 ++-- src/init.c | 16 +++++++--------- src/io.c | 11 +++++++---- src/main.c | 4 ++-- src/navier-stokes.c | 52 ++++++++++++++++++++++++++++++++-------------------- 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;i0){ - 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;i0 ? -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; i0 ? -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)]); + } } -- cgit v1.2.3-54-g00ecf