diff options
| author | Ian Jauslin <ian@jauslin.org> | 2022-05-25 22:32:03 -0400 |
|---|---|---|
| committer | Ian Jauslin <ian@jauslin.org> | 2022-05-25 22:32:03 -0400 |
| commit | fa8dc681cc502f3111f7e0c7bcab6fbe20013971 (patch) | |
| tree | 8ec6f64736f04688126b06f11dc1f6809ed23adb /src/navier-stokes.c | |
| parent | 5d92b56c01a79ecbb2df963b4732d0f12687cf75 (diff) | |
Enforce reality for whole convolution term
Diffstat (limited to 'src/navier-stokes.c')
| -rw-r--r-- | src/navier-stokes.c | 11 |
1 files changed, 4 insertions, 7 deletions
diff --git a/src/navier-stokes.c b/src/navier-stokes.c index e553022..fe2951f 100644 --- a/src/navier-stokes.c +++ b/src/navier-stokes.c @@ -100,7 +100,6 @@ int energy( ins_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3); if(t%print_freq==0){ - energy=0.; for(kx=-K1;kx<=K1;kx++){ for (ky=-K2;ky<=K2;ky++){ @@ -334,7 +333,6 @@ int ns_init_u( return 0; } - // next time step for Irreversible Navier-Stokes equation int ins_step( _Complex double* u, @@ -451,7 +449,8 @@ int ins_rhs( 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/L/L*nu*(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]+(*g)(kx,ky)+4*M_PI*M_PI/L/L/sqrt(kx*kx+ky*ky)*ifft.fft[klookup(kx,ky,N1,N2)]; + // enforce the reality of u by adding ifft.fft(k) and the conjugate of ifft.fft(-k) + out[klookup(kx,ky,2*K1+1,2*K2+1)]=-4*M_PI*M_PI/L/L*nu*(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]+(*g)(kx,ky)+4*M_PI*M_PI/L/L/sqrt(kx*kx+ky*ky)*(ifft.fft[klookup(kx,ky,N1,N2)]+conj(ifft.fft[klookup(-kx,-ky,N1,N2)]))/2; } } } @@ -495,8 +494,7 @@ int ns_T( fftw_execute(fft2.fft_plan); // write to ifft for(i=0;i<N1*N2;i++){ - // control numerical truncation by taking imaginary part (fft1 and fft2 should be purely imaginary) - ifft.fft[i]=-(__imag__ fft1.fft[i])*(__imag__ fft2.fft[i]); + ifft.fft[i]=fft1.fft[i]*fft2.fft[i]; } // F(py/|p|*u)*F(qx*|q|*u) @@ -520,8 +518,7 @@ int ns_T( fftw_execute(fft2.fft_plan); // write to ifft for(i=0;i<N1*N2;i++){ - // control numerical truncation by taking imaginary part (fft1 and fft2 should be purely imaginary) - ifft.fft[i]=ifft.fft[i]+(__imag__ fft1.fft[i])*(__imag__ fft2.fft[i]); + ifft.fft[i]=ifft.fft[i]-fft1.fft[i]*fft2.fft[i]; } // inverse fft |