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

Enforce T(u,-k)=T(u,k)^* in ns_T

Browse Source
This commit is contained in:
Ian Jauslin 2023-03-31 15:56:50 -04:00
parent d69f386547
commit 2958cc0313
1 changed files with 15 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
src/navier-stokes.c
View File

@ -383,8 +383,7 @@ 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){
// 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[klookup(kx,ky,2*K1+1,2*K2+1)]+4*M_PI*M_PI/L/L/sqrt(kx*kx+ky*ky)*ifft.fft[klookup(kx,ky,N1,N2)];
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[klookup(kx,ky,2*K1+1,2*K2+1)]+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;
} }
} }
} }
@ -458,6 +457,13 @@ int ns_T(
// inverse fft // inverse fft
fftw_execute(ifft.fft_plan); fftw_execute(ifft.fft_plan);
// enforce T(u,-k)=T(u,k)^*
for(kx=-K1;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){
ifft.fft[klookup(kx,ky,N1,N2)]=(ifft.fft[klookup(kx,ky,N1,N2)]+conj(ifft.fft[klookup(-kx,-ky,N1,N2)]))/2;
}
}
return(0); return(0);
} }
@ -498,6 +504,13 @@ int ns_T_nofft(
} }
} }
// enforce T(u,-k)=T(u,k)^*
for(kx=-K1;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){
out[klookup(kx,ky,N1,N2)]=(out[klookup(kx,ky,N1,N2)]+conj(out[klookup(-kx,-ky,N1,N2)]))/2;
}
}
return 0; return 0;
} }