diff options
| -rw-r--r-- | src/main.c | 7 | ||||
| -rw-r--r-- | src/navier-stokes.c | 83 | ||||
| -rw-r--r-- | src/navier-stokes.h | 7 |
3 files changed, 97 insertions, 0 deletions
@@ -20,6 +20,7 @@ int set_parameter(char* lhs, char* rhs, int* K1, int* K2, int* N1, int* N2, unsi #define COMMAND_UK 1 #define COMMAND_ENSTROPHY 2 +#define COMMAND_QUIET 3 #define DRIVING_TEST 1 @@ -69,6 +70,9 @@ int main ( else if(command==COMMAND_ENSTROPHY){ enstrophy(K1, K2, N1, N2, nsteps, nu, delta, g, print_freq); } + else if(command==COMMAND_QUIET){ + quiet(K1, K2, N1, N2, nsteps, nu, delta, g); + } else if(command==0){ fprintf(stderr, "error: no command specified\n"); print_usage(); @@ -145,6 +149,9 @@ int read_args( else if(strcmp(argv[i], "enstrophy")==0){ *command=COMMAND_ENSTROPHY; } + else if(strcmp(argv[i], "quiet")==0){ + *command=COMMAND_QUIET; + } else{ fprintf(stderr, "error: unrecognized command: '%s'\n",argv[i]); return(-1); diff --git a/src/navier-stokes.c b/src/navier-stokes.c index 13d9531..87f5aba 100644 --- a/src/navier-stokes.c +++ b/src/navier-stokes.c @@ -105,6 +105,38 @@ int enstrophy( return(0); } +// compute solution as a function of time, but do not print anything (useful for debugging) +int quiet( + int K1, + int K2, + int N1, + int N2, + unsigned int nsteps, + double nu, + double delta, + _Complex double (*g)(int,int) +){ + _Complex double* u; + _Complex double* tmp1; + _Complex double* tmp2; + _Complex double* tmp3; + unsigned int t; + fft_vect fft1; + fft_vect fft2; + fft_vect ifft; + + ns_init_tmps(&u, &tmp1, &tmp2, &tmp3, &fft1, &fft2, &ifft, K1, K2, N1, N2); + ns_init_u(u, K1, K2); + + // iterate + for(t=0;t<nsteps;t++){ + ins_step(u, K1, K2, N1, N2, nu, delta, g, fft1, fft2, ifft, tmp1, tmp2, tmp3); + } + + ns_free_tmps(u, tmp1, tmp2, tmp3, fft1, fft2, ifft); + return(0); +} + // initialize vectors for computation int ns_init_tmps( @@ -332,6 +364,18 @@ int ins_rhs( // compute convolution term ns_T(u,K1,K2,N1,N2,fft1,fft2,ifft); + + /* + // compare convolution term (store result in fft1.fft) + ns_T_nofft(fft1.fft, u, K1, K2, N1, N2); + double cmp=0.; + for(i=0;i<N1*N2;i++){ + cmp+=(ifft.fft[i]-fft1.fft[i])*(ifft.fft[i]-fft1.fft[i]); + } + printf("% .15e\n",sqrt(cmp)); + */ + + for(i=0; i<(2*K1+1)*(2*K2+1); i++){ out[i]=0; } @@ -415,6 +459,45 @@ int ns_T( return(0); } +// convolution term in right side of convolution equation, computed without fourier transform +int ns_T_nofft( + _Complex double* out, + _Complex double* u, + int K1, + int K2, + int N1, + int N2 +){ + int kx,ky; + int px,py; + int qx,qy; + + // loop over K's (needs N1>=4*K1+1 and N2>=4*K2+1) + if (N1<4*K1+1 || N2<4*K2+1){ + fprintf(stderr,"error: N1 and N2 need t be >= 4*K1+1 and 4*K2+1 respectively\n"); + return(-1); + } + for(kx=-2*K1;kx<=2*K1;kx++){ + for(ky=-2*K2;ky<=2*K2;ky++){ + // init + out[klookup(kx,ky,N1,N2)]=0.; + + for(px=-K1;px<=K1;px++){ + for(py=-K2;py<=K2;py++){ + qx=kx-px; + qy=ky-py; + + // cutoff in q + 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)]; + } + } + } + } + } + + return 0; +} // compute alpha _Complex double compute_alpha( diff --git a/src/navier-stokes.h b/src/navier-stokes.h index 2ff4e3f..d144fce 100644 --- a/src/navier-stokes.h +++ b/src/navier-stokes.h @@ -17,6 +17,10 @@ int uk( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double d // compute enstrophy int enstrophy( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, _Complex double (*g)(int,int), unsigned int print_freq); +// compute solution as a function of time, but do not print anything (useful for debugging) +int quiet( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, _Complex double (*g)(int,int)); + + // initialize vectors for computation int ns_init_tmps( _Complex double **u, _Complex double ** tmp1, _Complex double **tmp2, _Complex double **tmp3, fft_vect* fft1, fft_vect *fft2, fft_vect *ifft, int K1, int K2, int N1, int N2); // release vectors @@ -34,6 +38,9 @@ int ins_rhs( _Complex double* out, _Complex double* u, int K1, int K2, int N1, i // 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); +// convolution term in right side of equation +int ns_T_nofft( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2); + // compute alpha _Complex double compute_alpha( _Complex double* u, int K1, int K2, _Complex double (*g)(int,int)); |