software/Nstrophy
1
Fork 0
Code Issues Pull Requests Releases Activity

Only store u[kx,ky] with kx>=0

Browse Source
This commit is contained in:
Ian Jauslin 2023-04-11 18:45:45 -04:00
parent bca217e698
commit d16c42d9f5
6 changed files with 82 additions and 64 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
docs/nstrophy_doc.tex
View File

@ -110,6 +110,7 @@ Thus,
. .
\label{realT} \label{realT}
\end{equation} \end{equation}
In order to keep the computation as quick as possible, we only compute and store the values for $k_1\geqslant 0$.
\bigskip \bigskip
\point{\bf FFT}. We compute T using a fast Fourier transform, defined as \point{\bf FFT}. We compute T using a fast Fourier transform, defined as

11
src/driving.c
View File

@ -9,16 +9,13 @@ int g_test(
int K2 int K2
){ ){
int kx,ky; int kx,ky;
for(kx=-K1;kx<=K1;kx++){ for(kx=0;kx<=K1;kx++){
for (ky=-K2;ky<=K2;ky++){ for (ky=-K2;ky<=K2;ky++){
if(kx==2 && ky==-1){ 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 if(kx==-2 && ky==1){
g[klookup(kx,ky,2*K1+1,2*K2+1)]=0.5-sqrt(3)/2*I;
} }
else{ 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 K2
){ ){
int i; 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.; g[i]=0.;
} }

15
src/init.c
View File

@ -26,8 +26,7 @@ int init_random (
if (kx!=0 || ky>0){ if (kx!=0 || ky>0){
x=-0.5+((double) rand())/RAND_MAX; x=-0.5+((double) rand())/RAND_MAX;
y=-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_sym(kx,ky,K2)]=x+y*I;
u0[klookup(-kx,-ky,2*K1+1,2*K2+1)]=conj(u0[klookup(kx,ky,2*K1+1,2*K2+1)]);
} }
} }
} }
@ -39,9 +38,9 @@ int init_random (
// fix the enstrophy // fix the enstrophy
rescale=compute_enstrophy(u0, K1, K2, L); 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++){ 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; int kx,ky;
double rescale; double rescale;
for(kx=-K1;kx<=K1;kx++){ for(kx=0;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){ 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 // fix the enstrophy
rescale=compute_enstrophy(u0, K1, K2, L); 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++){ 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);
} }
} }

13
src/io.c
View File

@ -12,9 +12,9 @@ int write_u(_Complex double* u, int K1, int K2, FILE* file){
return 0; return 0;
} }
for(kx=-K1;kx<=K1;kx++){ for(kx=0;kx<=K1;kx++){
for (ky=-K2;ky<=K2;ky++){ 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 // errors
if(ret!=4){ 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{ else{
if(kx>K1 || kx<-K1 || ky>K2 || ky<-K2){ 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{ else{
// set u // set u
u[klookup(kx, ky, 2*K1+1, 2*K2+1)]=r+i*I; u[klookup_sym(kx, ky, K2)]=r+i*I;
} }
} }
} }

102
src/navier-stokes.c
View File

@ -58,10 +58,10 @@ int uk(
for(kx=-K1;kx<=K1;kx++){ for(kx=-K1;kx<=K1;kx++){
for (ky=-K2;ky<=K2;ky++){ for (ky=-K2;ky<=K2;ky++){
if (kx*kx+ky*ky<=1){ 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"); fprintf(stderr,"\n");
@ -204,12 +204,12 @@ int ns_init_tmps(
unsigned int nthreads unsigned int nthreads
){ ){
// velocity field // 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 // allocate tmp vectors for computation
*tmp1=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),(2*K1+1)*(2*K2+1)); *tmp2=calloc(sizeof(_Complex double),(K1+1)*(2*K2+1));
*tmp3=calloc(sizeof(_Complex double),(2*K1+1)*(2*K2+1)); *tmp3=calloc(sizeof(_Complex double),(K1+1)*(2*K2+1));
// init threads // init threads
fftw_init_threads(); fftw_init_threads();
@ -266,7 +266,7 @@ int copy_u(
){ ){
int i; 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]; u[i]=u0[i];
} }
@ -297,61 +297,54 @@ int ns_step(
// k1 // k1
ns_rhs(tmp1, u, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible); ns_rhs(tmp1, u, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
// add to output // add to output
for(kx=-K1;kx<=K1;kx++){ for(kx=0;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){ 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 // u+h*k1/2
for(kx=-K1;kx<=K1;kx++){ for(kx=0;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){ 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 // k2
ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible); ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
// add to output // add to output
for(kx=-K1;kx<=K1;kx++){ for(kx=0;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){ 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 // u+h*k2/2
for(kx=-K1;kx<=K1;kx++){ for(kx=0;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){ 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 // k3
ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible); ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
// add to output // add to output
for(kx=-K1;kx<=K1;kx++){ for(kx=0;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){ 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 // u+h*k3
for(kx=-K1;kx<=K1;kx++){ for(kx=0;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){ 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 // k4
ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible); ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
// add to output // add to output
for(kx=-K1;kx<=K1;kx++){ for(kx=0;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){ 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)]; u[klookup_sym(kx,ky,K2)]=tmp3[klookup_sym(kx,ky,K2)]+delta/6*tmp1[klookup_sym(kx,ky,K2)];
}
}
// 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;
} }
} }
@ -397,13 +390,13 @@ int ns_rhs(
printf("% .15e\n",sqrt(cmp)); 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; out[i]=0;
} }
for(kx=-K1;kx<=K1;kx++){ for(kx=0;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){ for(ky=-K2;ky<=K2;ky++){
if(kx!=0 || ky!=0){ 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(kx=-K1;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){ for(ky=-K2;ky<=K2;ky++){
if(kx!=0 || ky!=0){ 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; 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)*u[klookup(kx,ky,2*K1+1,2*K2+1)]/N2; 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(kx=-K1;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){ for(ky=-K2;ky<=K2;ky++){
if(kx!=0 || ky!=0){ 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; 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)*u[klookup(kx,ky,2*K1+1,2*K2+1)]/N2; 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 // inverse fft
fftw_execute(ifft.fft_plan); fftw_execute(ifft.fft_plan);
/*
// enforce T(u,-k)=T(u,k)^* // enforce T(u,-k)=T(u,k)^*
for(kx=-K1;kx<=K1;kx++){ for(kx=-K1;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){ 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; 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); return(0);
} }
@ -519,7 +510,7 @@ int ns_T_nofft(
// cutoff in q // cutoff in q
if(qx>=-K1 && qx<=K1 && qy>=-K2 && qy<=K2 && qx*qx+qy*qy>0 && px*px+py*py>0){ 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(kx=-K1;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){ 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)]); 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)*u[klookup(kx,ky,2*K1+1,2*K2+1)]*conj(u[klookup(kx,ky,2*K1+1,2*K2+1)]); 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(kx=-K1;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){ 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)); 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)*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)); 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.; double out=0.;
for(kx=-K1;kx<=K1;kx++){ for(kx=-K1;kx<=K1;kx++){
for (ky=-K2;ky<=K2;ky++){ 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; return out;
@ -614,7 +605,7 @@ double compute_enstrophy(
double out=0.; double out=0.;
for(kx=-K1;kx<=K1;kx++){ for(kx=-K1;kx<=K1;kx++){
for (ky=-K2;ky<=K2;ky++){ 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; return out;
@ -630,3 +621,26 @@ int klookup(
){ ){
return (kx>=0 ? kx : S1+kx)*S2 + (ky>=0 ? ky : S2+ky); 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)]));
}

4
src/navier-stokes.h
View File

@ -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 // get index for kx,ky in array of size S
int klookup( int kx, int ky, int S1, int S2); 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 #endif