diff options
Diffstat (limited to 'src/navier-stokes.c')
-rw-r--r-- | src/navier-stokes.c | 38 |
1 files changed, 16 insertions, 22 deletions
diff --git a/src/navier-stokes.c b/src/navier-stokes.c index 3dd7484..86a8497 100644 --- a/src/navier-stokes.c +++ b/src/navier-stokes.c @@ -68,7 +68,7 @@ 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 kx,ky; - // F(u/|p|)*F(q1*q2*u/|q|) + // F(px/|p|*u)*F(qy*|q|*u) // init to 0 for(kx=0; kx<params.N*params.N; kx++){ vects.fft1[kx]=0; @@ -79,11 +79,12 @@ int ins_rhs(_Complex double* out, _Complex double* u, ns_params params, fft_vect 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)]=u[KLOOKUP(kx,ky,params.S)]/sqrt(kx*kx+ky*ky); - vects.fft2[KLOOKUP(kx,ky,params.N)]=kx*ky*u[KLOOKUP(kx,ky,params.S)]/sqrt(kx*kx+ky*ky); + 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)]=ky*sqrt(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)]; } } } + // fft fftw_execute(vects.fft1_plan); fftw_execute(vects.fft2_plan); @@ -93,51 +94,44 @@ int ins_rhs(_Complex double* out, _Complex double* u, ns_params params, fft_vect 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); - - // 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*(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)]+params.g[KLOOKUP(kx,ky,params.S)]+4*M_PI*M_PI/sqrt(kx*kx+ky*ky)*vects.invfft[KLOOKUP(kx,ky,params.N)]*(kx*kx-ky*ky)/params.N/params.N; - } - } - } - // F(u/|p|)*F((q1*q1-q2*q2)*u/|q|) + // F(py/|p|*u)*F(qx*|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.fft2[KLOOKUP(kx,ky,params.N)]=(kx*kx-ky*ky)*u[KLOOKUP(kx,ky,params.S)]/sqrt(kx*kx+ky*ky); + vects.fft1[KLOOKUP(kx,ky,params.N)]=ky/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.invfft[KLOOKUP(kx,ky,params.N)]=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); // 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)]=out[KLOOKUP(kx,ky,params.S)]-4*M_PI*M_PI/sqrt(kx*kx+ky*ky)*vects.invfft[KLOOKUP(kx,ky,params.N)]*(kx*ky)/params.N/params.N; + 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)]+4*M_PI*M_PI/sqrt(kx*kx+ky*ky)*vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N; } } } |