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

Do not enforce symmetry on T: only the k>0 matter anyways

Browse Source
This commit is contained in:
Ian Jauslin
2023-04-22 15:01:31 -04:00
parent 731244a83d
commit 5d0a1bcc6f
3 changed files with 9 additions and 22 deletions
Download Patch File Download Diff File Expand all files Collapse all files

27
src/navier-stokes.c
View File

@@ -207,6 +207,7 @@ int quiet(
double nu,
double delta,
double L,
uint64_t starting_time,
_Complex double* u0,
_Complex double* g,
bool irreversible,
@@ -227,7 +228,7 @@ int quiet(
copy_u(u, u0, K1, K2);
// iterate
for(t=0;nsteps==0 || t<nsteps;t++){
for(t=starting_time;nsteps==0 || t<starting_time+nsteps;t++){
ns_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
}
@@ -511,13 +512,6 @@ int ns_T(
// inverse fft
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);
}
@@ -534,13 +528,13 @@ int ns_T_nofft(
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");
// loop over K's (needs N1>=2*K1+1 and N2>=2*K2+1)
if (N1<2*K1+1 || N2<2*K2+1){
fprintf(stderr,"error: N1 and N2 need t be >= 2*K1+1 and 2*K2+1 respectively\n");
return(-1);
}
for(kx=-2*K1;kx<=2*K1;kx++){
for(ky=-2*K2;ky<=2*K2;ky++){
for(kx=-K1;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){
// init
out[klookup(kx,ky,N1,N2)]=0.;
@@ -558,13 +552,6 @@ 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;
}