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

Enforce reality for whole convolution term

Browse Source
This commit is contained in:
Ian Jauslin 2022-05-25 22:32:03 -04:00
parent 5d92b56c01
commit fa8dc681cc
2 changed files with 5 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/main.c
View File

@ -77,7 +77,7 @@ int main (
energy(K1, K2, N1, N2, nsteps, nu, delta, L, g, print_freq, nthreads);
}
else if(command==COMMAND_ENSTROPHY){
energy(K1, K2, N1, N2, nsteps, nu, delta, L, g, print_freq, nthreads);
enstrophy(K1, K2, N1, N2, nsteps, nu, delta, L, g, print_freq, nthreads);
}
else if(command==COMMAND_QUIET){
quiet(K1, K2, N1, N2, nsteps, nu, delta, L, g, nthreads);

11
src/navier-stokes.c
View File

@ -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