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

compute convolution term in its own function

Browse Source
This commit is contained in:
Ian Jauslin 2022-05-18 10:53:00 +02:00
parent 232822a183
commit eca5070274
2 changed files with 38 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

51
src/navier-stokes.c
View File

@ -329,6 +329,37 @@ int ins_rhs(
int kx,ky; int kx,ky;
int i; int i;
// compute convolution term
ns_T(u,K1,K2,N1,N2,fft1,fft2,ifft);
for(i=0; i<(2*K1+1)*(2*K2+1); i++){
out[i]=0;
}
for(kx=-K1;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*nu*(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]+(*g)(kx,ky)+4*M_PI*M_PI/sqrt(kx*kx+ky*ky)*ifft.fft[klookup(kx,ky,N1,N2)];
}
}
}
return(0);
}
// convolution term in right side of convolution equation
int ns_T(
_Complex double* u,
int K1,
int K2,
int N1,
int N2,
fft_vect fft1,
fft_vect fft2,
fft_vect ifft
){
int kx,ky;
int i;
// F(px/|p|*u)*F(qy*|q|*u) // F(px/|p|*u)*F(qy*|q|*u)
// init to 0 // init to 0
for(i=0; i<N1*N2; i++){ for(i=0; i<N1*N2; i++){
@ -340,8 +371,8 @@ int ins_rhs(
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)]; 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)]; fft2.fft[klookup(kx,ky,N1,N2)]=ky*sqrt(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]/N2;
} }
} }
} }
@ -364,8 +395,8 @@ int ins_rhs(
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)]; 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)]; fft2.fft[klookup(kx,ky,N1,N2)]=kx*sqrt(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]/N2;
} }
} }
} }
@ -381,18 +412,6 @@ int ins_rhs(
// inverse fft // inverse fft
fftw_execute(ifft.fft_plan); fftw_execute(ifft.fft_plan);
// write out
for(i=0; i<(2*K1+1)*(2*K2+1); i++){
out[i]=0;
}
for(kx=-K1;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*nu*(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]+(*g)(kx,ky)+4*M_PI*M_PI/sqrt(kx*kx+ky*ky)*ifft.fft[klookup(kx,ky,N1,N2)]/N1/N2;
}
}
}
return(0); return(0);
} }

3
src/navier-stokes.h
View File

@ -31,6 +31,9 @@ int ins_step( _Complex double* u, int K1, int K2, int N1, int N2, double nu, dou
// right side of Irreversible Navier-Stokes equation // right side of Irreversible Navier-Stokes equation
int ins_rhs( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2, double nu, _Complex double (*g)(int,int), fft_vect fft1, fft_vect fft2, fft_vect ifft); int ins_rhs( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2, double nu, _Complex double (*g)(int,int), fft_vect fft1, fft_vect fft2, fft_vect ifft);
// convolution term in right side of equation
int ns_T( _Complex double* u, int K1, int K2, int N1, int N2, fft_vect fft1, fft_vect fft2, fft_vect ifft);
// compute alpha // compute alpha
_Complex double compute_alpha( _Complex double* u, int K1, int K2, _Complex double (*g)(int,int)); _Complex double compute_alpha( _Complex double* u, int K1, int K2, _Complex double (*g)(int,int));