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

Remove old ways of computing convolution

Browse Source
This commit is contained in:
Ian Jauslin 2018-01-11 22:48:55 +00:00
parent 01f47ace67
commit 7ee5507b93
1 changed files with 0 additions and 136 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

136
src/navier-stokes.c
View File

@ -3,8 +3,6 @@
#define M_PI 3.14159265358979323846 #define M_PI 3.14159265358979323846
#define CHK 1
// next time step for Irreversible Navier-Stokes equation // next time step for Irreversible Navier-Stokes equation
int ins_step(_Complex double* u, ns_params params, fft_vects vects, _Complex double* tmp1, _Complex double* tmp2, _Complex double* tmp3){ int ins_step(_Complex double* u, ns_params params, fft_vects vects, _Complex double* tmp1, _Complex double* tmp2, _Complex double* tmp3){
int kx,ky; int kx,ky;
@ -70,7 +68,6 @@ int ins_step(_Complex double* u, ns_params params, fft_vects vects, _Complex dou
int ins_rhs(_Complex double* out, _Complex double* u, ns_params params, fft_vects vects){ int ins_rhs(_Complex double* out, _Complex double* u, ns_params params, fft_vects vects){
int kx,ky; int kx,ky;
#if CHK==0
// F(u/|p|)*F(q1*q2*u/|q|) // F(u/|p|)*F(q1*q2*u/|q|)
// init to 0 // init to 0
for(kx=0; kx<params.N*params.N; kx++){ for(kx=0; kx<params.N*params.N; kx++){
@ -144,140 +141,7 @@ int ins_rhs(_Complex double* out, _Complex double* u, ns_params params, fft_vect
} }
} }
} }
#elif CHK == 1
// F(-p2/|p|*u)*F(q1*|q|*u)
// init to 0
for(kx=0; kx<params.N*params.N; kx++){
vects.fft1[kx]=0;
vects.fft2[kx]=0;
vects.invfft[kx]=0;
}
// fill modes
for(kx=-params.K;kx<=params.K;kx++){
for(ky=-params.K;ky<=params.K;ky++){
if(kx!=0 || ky!=0){
vects.fft1[KLOOKUP(kx,ky,params.N)]=-kx/sqrt(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)];
vects.fft2[KLOOKUP(kx,ky,params.N)]=kx*sqrt(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)];
}
}
}
// fft
fftw_execute(vects.fft1_plan);
fftw_execute(vects.fft2_plan);
// write to invfft
for(kx=-2*params.K;kx<=2*params.K;kx++){
for(ky=-2*params.K;ky<=2*params.K;ky++){
vects.invfft[KLOOKUP(kx,ky,params.N)]=vects.fft1[KLOOKUP(kx,ky,params.N)]*vects.fft2[KLOOKUP(kx,ky,params.N)] - vects.fft1[KLOOKUP(ky,kx,params.N)]*vects.fft2[KLOOKUP(ky,kx,params.N)];
}
}
// inverse fft
fftw_execute(vects.invfft_plan);
// write out
for(kx=0; kx<params.S*params.S; kx++){
out[kx]=0;
}
for(kx=-params.K;kx<=params.K;kx++){
for(ky=-params.K;ky<=params.K;ky++){
if(kx!=0 || ky!=0){
out[KLOOKUP(kx,ky,params.S)]=-4*M_PI*M_PI*params.nu/params.S*(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)]+params.g[KLOOKUP(kx,ky,params.S)]+2*M_PI/sqrt(kx*kx+ky*ky)/params.S*vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N;
}
}
}
#elif CHK==2
// F(u)*F(q1*u)
// init to 0
for(kx=0; kx<params.N*params.N; kx++){
vects.fft1[kx]=0;
vects.fft2[kx]=0;
vects.invfft[kx]=0;
}
// fill modes
for(kx=-params.K;kx<=params.K;kx++){
for(ky=-params.K;ky<=params.K;ky++){
// u_k
vects.fft1[KLOOKUP(kx,ky,params.N)]=u[KLOOKUP(kx,ky,params.S)];
// 2i\pi k_x u_k
vects.fft2[KLOOKUP(kx,ky,params.N)]=2*M_PI*I*kx*u[KLOOKUP(kx,ky,params.S)];
}
}
// fft
fftw_execute(vects.fft1_plan);
fftw_execute(vects.fft2_plan);
// write to invfft
for(kx=-2*params.K;kx<=2*params.K;kx++){
for(ky=-2*params.K;ky<=2*params.K;ky++){
vects.invfft[KLOOKUP(kx,ky,params.N)]=vects.fft1[KLOOKUP(kx,ky,params.N)]*vects.fft2[KLOOKUP(kx,ky,params.N)];
}
}
// F(p1/p2*u)*F(q2*u)
// init to 0
for(kx=0; kx<params.N*params.N; kx++){
vects.fft1[kx]=0;
vects.fft2[kx]=0;
}
// fill modes
for(kx=-params.K;kx<=params.K;kx++){
for(ky=-params.K;ky<=params.K;ky++){
// k_x/k_y u_k
if(ky!=0){
vects.fft1[KLOOKUP(kx,ky,params.N)]=kx/ky*u[KLOOKUP(kx,ky,params.S)];
}
// 2i\pi k_y u_k
vects.fft2[KLOOKUP(kx,ky,params.N)]=2*M_PI*I*ky*u[KLOOKUP(kx,ky,params.S)];
}
}
// fft
fftw_execute(vects.fft1_plan);
fftw_execute(vects.fft2_plan);
// write to invfft
for(kx=-2*params.K;kx<=2*params.K;kx++){
for(ky=-2*params.K;ky<=2*params.K;ky++){
vects.invfft[KLOOKUP(kx,ky,params.N)]+=-vects.fft1[KLOOKUP(kx,ky,params.N)]*vects.fft2[KLOOKUP(kx,ky,params.N)];
}
}
// inverse fft
fftw_execute(vects.invfft_plan);
/*
// check: convolution instead of fft
for(kx=0; kx<params.S*params.S; kx++){
out[kx]=0;
}
for(kx=-params.K;kx<=params.K;kx++){
for(ky=-params.K;ky<=params.K;ky++){
for(px=-params.K;px<=params.K;px++){
for(py=-params.K;py<=params.K;py++){
if(kx-px<=params.K && kx-px>=-params.K && ky-py<=params.K && ky-py>=-params.K){
out[KLOOKUP(kx,ky,params.S)]+=2*M_PI*I*(u[KLOOKUP(px,py,params.S)]*(kx-px)*u[KLOOKUP(kx-px,ky-py,params.S)]-(py==0?0:px/py*u[KLOOKUP(px,py,params.S)]*(ky-py)*u[KLOOKUP(kx-px,ky-py,params.S)]));
}
}
}
dd=(__real__ vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N-__real__ out[KLOOKUP(kx,ky,params.S)])*(__real__ vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N-__real__ out[KLOOKUP(kx,ky,params.S)])+(__imag__ vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N-__imag__ out[KLOOKUP(kx,ky,params.S)])*(__imag__ vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N-__imag__ out[KLOOKUP(kx,ky,params.S)]);
if(dd>1e-25){
printf("%d %d % .8e\n",kx,ky, dd);
}
}
}
*/
// write out
for(kx=0; kx<params.S*params.S; kx++){
out[kx]=0;
}
for(kx=-params.K;kx<=params.K;kx++){
for(ky=-params.K;ky<=params.K;ky++){
out[KLOOKUP(kx,ky,params.S)]=-4*M_PI*M_PI*params.nu*(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)]+params.g[KLOOKUP(kx,ky,params.S)]+vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N;
}
}
#endif
return(0); return(0);
} }