From d16c42d9f5a40b94406a859fa510bba96480d5e8 Mon Sep 17 00:00:00 2001 From: Ian Jauslin Date: Tue, 11 Apr 2023 18:45:45 -0400 Subject: Only store u[kx,ky] with kx>=0 --- src/driving.c | 11 +++--- src/init.c | 15 ++++---- src/io.c | 13 ++++--- src/navier-stokes.c | 102 +++++++++++++++++++++++++++++----------------------- src/navier-stokes.h | 4 +++ 5 files changed, 81 insertions(+), 64 deletions(-) (limited to 'src') 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 -- cgit v1.2.3-54-g00ecf