Remove term from alpha that is equal to 0 anyways

This commit is contained in:
Ian Jauslin 2023-04-21 12:25:45 -04:00
parent 063dae642a
commit bc05a28b30
2 changed files with 5 additions and 30 deletions

View File

@ -132,7 +132,7 @@ int eea(
ns_T(u,K1,K2,N1,N2,fft1,fft2,ifft); 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(u, K1, K2, g, L);
enstrophy=compute_enstrophy(u, K1, K2, L); enstrophy=compute_enstrophy(u, K1, K2, L);
// running average // running average
@ -431,7 +431,7 @@ int ns_rhs(
if (irreversible) { if (irreversible) {
alpha=nu; alpha=nu;
} else { } else {
alpha=compute_alpha(u,K1,K2,N1,N2,g,L,ifft.fft); alpha=compute_alpha(u,K1,K2,g,L);
} }
/* /*
@ -586,11 +586,8 @@ double compute_alpha(
_Complex double* u, _Complex double* u,
int K1, int K1,
int K2, int K2,
int N1,
int N2,
_Complex double* g, _Complex double* g,
double L, double L
_Complex double* T
){ ){
_Complex double num=0; _Complex double num=0;
double denom=0; double denom=0;
@ -601,7 +598,7 @@ double compute_alpha(
for(kx=-K1;kx<=K1;kx++){ for(kx=-K1;kx<=K1;kx++){
for(ky=-K2;ky<=K2;ky++){ for(ky=-K2;ky<=K2;ky++){
num+=(L*L/4/M_PI/M_PI*(kx*kx+ky*ky)*getval_sym(g, kx,ky,K2)+sqrt(kx*kx+ky*ky)*T[klookup(kx,ky,N1,N2)])*conj(getval_sym(u, kx,ky,K2)); 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));
denom+=__real__ (kx*kx+ky*ky)*(kx*kx+ky*ky)*getval_sym(u, kx,ky,K2)*conj(getval_sym(u, kx,ky,K2)); denom+=__real__ (kx*kx+ky*ky)*(kx*kx+ky*ky)*getval_sym(u, kx,ky,K2)*conj(getval_sym(u, kx,ky,K2));
} }
} }
@ -610,28 +607,6 @@ double compute_alpha(
} }
/*
// compute alpha
double compute_alpha(
_Complex double* u,
int K1,
int K2,
_Complex double* g
){
_Complex double num=0;
_Complex double denom=0;
int kx,ky;
for(kx=-K1;kx<=K1;kx++){
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+=(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));
}
}
return __real__ (num/denom);
}*/
// compute energy // compute energy
double compute_energy( double compute_energy(
_Complex double* u, _Complex double* u,

View File

@ -48,7 +48,7 @@ int ns_T( _Complex double* u, int K1, int K2, int N1, int N2, fft_vect fft1, fft
int ns_T_nofft( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2); int ns_T_nofft( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2);
// 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, _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