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/navier-stokes.c | 102 +++++++++++++++++++++++++++++----------------------- 1 file changed, 58 insertions(+), 44 deletions(-) (limited to 'src/navier-stokes.c') 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)])); +} -- cgit v1.2.3-54-g00ecf