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

Quiet mode

Browse Source
This commit is contained in:
Ian Jauslin 2022-05-18 22:22:42 +02:00
parent eca5070274
commit 199b8f0df5
3 changed files with 97 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
src/main.c
View File

@ -20,6 +20,7 @@ int set_parameter(char* lhs, char* rhs, int* K1, int* K2, int* N1, int* N2, unsi
#define COMMAND_UK 1
#define COMMAND_ENSTROPHY 2
#define COMMAND_QUIET 3
#define DRIVING_TEST 1
@ -69,6 +70,9 @@ int main (
else if(command==COMMAND_ENSTROPHY){
enstrophy(K1, K2, N1, N2, nsteps, nu, delta, g, print_freq);
}
else if(command==COMMAND_QUIET){
quiet(K1, K2, N1, N2, nsteps, nu, delta, g);
}
else if(command==0){
fprintf(stderr, "error: no command specified\n");
print_usage();
@ -145,6 +149,9 @@ int read_args(
else if(strcmp(argv[i], "enstrophy")==0){
*command=COMMAND_ENSTROPHY;
}
else if(strcmp(argv[i], "quiet")==0){
*command=COMMAND_QUIET;
}
else{
fprintf(stderr, "error: unrecognized command: '%s'\n",argv[i]);
return(-1);

83
src/navier-stokes.c
View File

@ -105,6 +105,38 @@ int enstrophy(
return(0);
}
// compute solution as a function of time, but do not print anything (useful for debugging)
int quiet(
int K1,
int K2,
int N1,
int N2,
unsigned int nsteps,
double nu,
double delta,
_Complex double (*g)(int,int)
){
_Complex double* u;
_Complex double* tmp1;
_Complex double* tmp2;
_Complex double* tmp3;
unsigned int t;
fft_vect fft1;
fft_vect fft2;
fft_vect ifft;
ns_init_tmps(&u, &tmp1, &tmp2, &tmp3, &fft1, &fft2, &ifft, K1, K2, N1, N2);
ns_init_u(u, K1, K2);
// iterate
for(t=0;t<nsteps;t++){
ins_step(u, K1, K2, N1, N2, nu, delta, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
}
ns_free_tmps(u, tmp1, tmp2, tmp3, fft1, fft2, ifft);
return(0);
}
// initialize vectors for computation
int ns_init_tmps(
@ -332,6 +364,18 @@ int ins_rhs(
// compute convolution term
ns_T(u,K1,K2,N1,N2,fft1,fft2,ifft);
/*
// 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<(2*K1+1)*(2*K2+1); i++){
out[i]=0;
}
@ -415,6 +459,45 @@ int ns_T(
return(0);
}
// convolution term in right side of convolution equation, computed without fourier transform
int ns_T_nofft(
_Complex double* out,
_Complex double* u,
int K1,
int K2,
int N1,
int N2
){
int kx,ky;
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");
return(-1);
}
for(kx=-2*K1;kx<=2*K1;kx++){
for(ky=-2*K2;ky<=2*K2;ky++){
// init
out[klookup(kx,ky,N1,N2)]=0.;
for(px=-K1;px<=K1;px++){
for(py=-K2;py<=K2;py++){
qx=kx-px;
qy=ky-py;
// cutoff in q
if(qx>=-K1 && qx<=K1 && qy>=-K2 && qy<=K2 && qx*qx+qy*qy>0 && px*px+py*py>0){
out[klookup(kx,ky,N1,N2)]+=(-qx*py+qy*px)*sqrt(qx*qx+qy*qy)/sqrt(px*px+py*py)*u[klookup(px,py,2*K1+1,2*K2+1)]*u[klookup(qx,qy,2*K1+1,2*K2+1)];
}
}
}
}
}
return 0;
}
// compute alpha
_Complex double compute_alpha(

7
src/navier-stokes.h
View File

@ -17,6 +17,10 @@ int uk( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double d
// compute enstrophy
int enstrophy( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, _Complex double (*g)(int,int), unsigned int print_freq);
// compute solution as a function of time, but do not print anything (useful for debugging)
int quiet( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, _Complex double (*g)(int,int));
// initialize vectors for computation
int ns_init_tmps( _Complex double **u, _Complex double ** tmp1, _Complex double **tmp2, _Complex double **tmp3, fft_vect* fft1, fft_vect *fft2, fft_vect *ifft, int K1, int K2, int N1, int N2);
// release vectors
@ -34,6 +38,9 @@ int ins_rhs( _Complex double* out, _Complex double* u, int K1, int K2, int N1, i
// convolution term in right side of equation
int ns_T( _Complex double* u, int K1, int K2, int N1, int N2, fft_vect fft1, fft_vect fft2, fft_vect ifft);
// convolution term in right side of equation
int ns_T_nofft( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2);
// compute alpha
_Complex double compute_alpha( _Complex double* u, int K1, int K2, _Complex double (*g)(int,int));