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

Only use T term in alpha for ns_step computation

Browse Source
This commit is contained in:
Ian Jauslin 2023-04-21 17:36:17 -04:00
parent a53aec9b22
commit 731244a83d
2 changed files with 16 additions and 23 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

37
src/navier-stokes.c
View File

@ -128,11 +128,8 @@ int eea(
for(t=starting_time;nsteps==0 || t<starting_time+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); ns_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
// convolution term
ns_T(u,K1,K2,N1,N2,fft1,fft2,ifft);
energy=compute_energy(u, K1, K2); energy=compute_energy(u, K1, K2);
alpha=compute_alpha(u, K1, K2, N1, N2, g, L, ifft.fft); alpha=compute_alpha_fast(u, K1, K2, g, L);
enstrophy=compute_enstrophy(u, K1, K2, L); enstrophy=compute_enstrophy(u, K1, K2, L);
// running average // running average
@ -434,16 +431,6 @@ int ns_rhs(
alpha=compute_alpha(u,K1,K2,N1,N2,g,L,ifft.fft); alpha=compute_alpha(u,K1,K2,N1,N2,g,L,ifft.fft);
} }
/*
// 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<(K1+1)*(2*K2+1); i++){ for(i=0; i<(K1+1)*(2*K2+1); i++){
out[i]=0; out[i]=0;
} }
@ -609,28 +596,32 @@ double compute_alpha(
return __real__ num/denom; return __real__ num/denom;
} }
/*
// compute alpha // compute alpha
double compute_alpha( // do not include the term involving T, which, in the limit K->infty, vanishes
double compute_alpha_fast(
_Complex double* u, _Complex double* u,
int K1, int K1,
int K2, int K2,
_Complex double* g _Complex double* g,
double L
){ ){
_Complex double num=0; _Complex double num=0;
_Complex double denom=0; double denom=0;
int kx,ky; int kx,ky;
num=0.;
denom=0.;
for(kx=-K1;kx<=K1;kx++){ for(kx=-K1;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){ for(ky=-K2;ky<=K2;ky++){
denom+=(kx*kx+ky*ky)*(kx*kx+ky*ky)*getval_sym(u, kx,ky,K2)*conj(getval_sym(u, kx,ky,K2))*(1+(ky!=0?kx*kx/ky/ky:0)); num+=L*L/4/M_PI/M_PI*(kx*kx+ky*ky)*getval_sym(g, kx,ky,K2)*conj(getval_sym(u, kx,ky,K2));
num+=(kx*kx+ky*ky)*getval_sym(u, kx,ky,K2)*conj(g[getval(kx,ky,K1,K2)])*(1+(ky!=0?kx*kx/ky/ky:0)); denom+=__real__ (kx*kx+ky*ky)*(kx*kx+ky*ky)*getval_sym(u, kx,ky,K2)*conj(getval_sym(u, kx,ky,K2));
} }
} }
return __real__ (num/denom); return __real__ num/denom;
}*/ }
// compute energy // compute energy
double compute_energy( double compute_energy(

2
src/navier-stokes.h
View File

@ -49,6 +49,8 @@ int ns_T_nofft( _Complex double* out, _Complex double* u, int K1, int K2, int N1
// compute alpha // compute alpha
double compute_alpha( _Complex double* u, int K1, int K2, int N1, int N2, _Complex double* g, double L, _Complex double* T); double compute_alpha( _Complex double* u, int K1, int K2, int N1, int N2, _Complex double* g, double L, _Complex double* T);
// do not include the term involving T, which, in the limit K->infty, vanishes
double compute_alpha_fast( _Complex double* u, int K1, int K2, _Complex double* g, double L);
// compute energy // compute energy
double compute_energy( _Complex double* u, int K1, int K2); double compute_energy( _Complex double* u, int K1, int K2);
// compute enstrophy // compute enstrophy