diff options
author | Ian Jauslin <ian@jauslin.org> | 2022-05-18 22:22:42 +0200 |
---|---|---|
committer | Ian Jauslin <ian@jauslin.org> | 2022-05-18 22:22:42 +0200 |
commit | 199b8f0df5adeaaac4ca2afc2eff5237dfee36c3 (patch) | |
tree | dfcd4ba2d05757070fe2883b8c64c84a8dbd67a4 /src/navier-stokes.c | |
parent | eca50702746fc0e8c933bac32cd4e5623d88ca53 (diff) |
Quiet mode
Diffstat (limited to 'src/navier-stokes.c')
-rw-r--r-- | src/navier-stokes.c | 83 |
1 files changed, 83 insertions, 0 deletions
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( |