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

Remove terms with kx=0 and ky<=0

Browse Source
This commit is contained in:
Ian Jauslin 2023-04-25 17:51:14 -04:00
parent 519bb72726
commit d838060265
5 changed files with 50 additions and 37 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/driving.c
View File

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

10
src/init.c
View File

@ -22,14 +22,12 @@ int init_random (
// random init (set half, then the other half are the conjugates) // random init (set half, then the other half are the conjugates)
for(kx=0;kx<=K1;kx++){ for(kx=0;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){ for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
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_sym(kx,ky,K2)]=x+y*I; u0[klookup_sym(kx,ky,K2)]=x+y*I;
} }
} }
}
if (irreversible) { if (irreversible) {
// fix the energy // fix the energy
@ -39,7 +37,7 @@ int init_random (
rescale=compute_enstrophy(u0, K1, K2, L); rescale=compute_enstrophy(u0, K1, K2, L);
} }
for(kx=0;kx<=K1;kx++){ 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); 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; double rescale;
for(kx=0;kx<=K1;kx++){ 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)); 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); rescale=compute_enstrophy(u0, K1, K2, L);
} }
for(kx=0;kx<=K1;kx++){ 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); u0[klookup_sym(kx,ky,K2)]=u0[klookup_sym(kx,ky,K2)]*sqrt(init_en/rescale);
} }
} }

11
src/io.c
View File

@ -14,7 +14,7 @@ int write_vec(_Complex double* vec, int K1, int K2, FILE* file){
} }
for(kx=0;kx<=K1;kx++){ 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)]); 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; 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; 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); fprintf(stderr, "warning: reading line %d: kx or ky out of bounds: %d,%d, skipping\n", counter, kx, ky);
} }
else if (kx<0){ 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{ else{
// set output // 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; return 0;
} }

4
src/main.c
View File

@ -635,7 +635,7 @@ int set_parameter(
_Complex double* set_driving( _Complex double* set_driving(
nstrophy_parameters parameters 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){ switch(parameters.driving){
case DRIVING_ZERO: case DRIVING_ZERO:
@ -664,7 +664,7 @@ _Complex double* set_driving(
_Complex double* set_init( _Complex double* set_init(
nstrophy_parameters parameters 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){ switch(parameters.init){
case INIT_RANDOM: case INIT_RANDOM:

50
src/navier-stokes.c
View File

@ -241,12 +241,12 @@ int ns_init_tmps(
unsigned int nthreads unsigned int nthreads
){ ){
// velocity field // 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 // allocate tmp vectors for computation
*tmp1=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+1)*(2*K2+1)); *tmp2=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
*tmp3=calloc(sizeof(_Complex double),(K1+1)*(2*K2+1)); *tmp3=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
// init threads // init threads
fftw_init_threads(); fftw_init_threads();
@ -303,7 +303,7 @@ int copy_u(
){ ){
int i; int i;
for(i=0;i<(K1+1)*(2*K2+1);i++){ for(i=0;i<K1*(2*K2+1)+K2;i++){
u[i]=u0[i]; u[i]=u0[i];
} }
@ -335,14 +335,14 @@ int ns_step(
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=0;kx<=K1;kx++){ 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)]=u[klookup_sym(kx,ky,K2)]+delta/6*tmp1[klookup_sym(kx,ky,K2)]; 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=0;kx<=K1;kx++){ 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)]; 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); ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
// add to output // add to output
for(kx=0;kx<=K1;kx++){ 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)]; tmp3[klookup_sym(kx,ky,K2)]+=delta/3*tmp1[klookup_sym(kx,ky,K2)];
} }
} }
// u+h*k2/2 // u+h*k2/2
for(kx=0;kx<=K1;kx++){ 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)]; 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); ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
// add to output // add to output
for(kx=0;kx<=K1;kx++){ 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)]; tmp3[klookup_sym(kx,ky,K2)]+=delta/3*tmp1[klookup_sym(kx,ky,K2)];
} }
} }
// u+h*k3 // u+h*k3
for(kx=0;kx<=K1;kx++){ 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)]; 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); ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
// add to output // add to output
for(kx=0;kx<=K1;kx++){ 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)]; u[klookup_sym(kx,ky,K2)]=tmp3[klookup_sym(kx,ky,K2)]+delta/6*tmp1[klookup_sym(kx,ky,K2)];
} }
} }
@ -417,16 +417,14 @@ int ns_rhs(
alpha=compute_alpha(u,K1,K2,g,L); alpha=compute_alpha(u,K1,K2,g,L);
} }
for(i=0; i<(K1+1)*(2*K2+1); i++){ for(i=0; i<K1*(2*K2+1)+K2; i++){
out[i]=0; out[i]=0;
} }
for(kx=0;kx<=K1;kx++){ for(kx=0;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){ for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
if(kx!=0 || ky!=0){
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)]; 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)];
} }
} }
}
return(0); return(0);
} }
@ -610,7 +608,7 @@ 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 // 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 klookup_sym(
int kx, int kx,
int ky, 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"); fprintf(stderr, "bug!: attempting to access a symmetrized vector at kx<0\n Contact Ian at ian.jauslin@rutgers.edu!\n");
exit(-1); 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^* // 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 ky,
int K2 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)]);
}
} }