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

Running average

Browse Source
This commit is contained in:
Ian Jauslin
2023-04-12 14:37:02 -04:00
parent d16c42d9f5
commit 59100a471f
3 changed files with 129 additions and 11 deletions
Download Patch File Download Diff File Expand all files Collapse all files

112
src/navier-stokes.c
View File

@@ -90,6 +90,7 @@ int eea(
_Complex double* g,
bool irreversible,
unsigned int print_freq,
unsigned int running_avg_window,
unsigned int starting_time,
unsigned int nthreads,
FILE* savefile
@@ -100,6 +101,11 @@ int eea(
_Complex double* tmp3;
double alpha, energy, enstrophy;
double avg_e,avg_a,avg_en;
// quantities needed to compute averages
double *save_print_e, *save_print_a, *save_print_en;
double *save_print_short_e, *save_print_short_a, *save_print_short_en;
// index
unsigned int i;
unsigned int t;
fft_vect fft1;
fft_vect fft2;
@@ -115,6 +121,28 @@ int eea(
avg_a=0;
avg_en=0;
// need this for running average
unsigned int short_len=running_avg_window % print_freq;
// number of full print intervals in running average window
unsigned int nr_print_in_window = (running_avg_window-short_len)/print_freq;
double avg_print_e=0;
double avg_print_a=0;
double avg_print_en=0;
double avg_print_short_e=0;
double avg_print_short_a=0;
double avg_print_short_en=0;
// save e a en for running average
// only need to save if window is non zero
if(running_avg_window!=0){
save_print_e=calloc(sizeof(double), nr_print_in_window);
save_print_a=calloc(sizeof(double), nr_print_in_window);
save_print_en=calloc(sizeof(double), nr_print_in_window);
save_print_short_e=calloc(sizeof(double), nr_print_in_window+1);
save_print_short_a=calloc(sizeof(double), nr_print_in_window+1);
save_print_short_en=calloc(sizeof(double), nr_print_in_window+1);
}
// iterate
for(t=starting_time;t<starting_time+nsteps;t++){
ns_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
@@ -127,18 +155,92 @@ int eea(
enstrophy=compute_enstrophy(u, K1, K2, L);
// running average
if(t>0){
avg_e=avg_e-(avg_e-energy)/t;
avg_a=avg_a-(avg_a-alpha)/t;
avg_en=avg_en-(avg_en-enstrophy)/t;
// if window=0, then take average over all times
if(running_avg_window==0){
if(t!=0){
avg_e=avg_e-(avg_e-energy)/t;
avg_a=avg_a-(avg_a-alpha)/t;
avg_en=avg_en-(avg_en-enstrophy)/t;
}
}
else{
// reset print interval averages
if(t % print_freq == 1){
avg_print_e=0;
avg_print_a=0;
avg_print_en=0;
avg_print_short_e=0;
avg_print_short_a=0;
avg_print_short_en=0;
}
// compute average over print interval
avg_print_e=avg_print_e+energy/print_freq;
avg_print_a=avg_print_a+alpha/print_freq;
avg_print_en=avg_print_en+enstrophy/print_freq;
// compute average over the last short_len elements of a print interval
if (short_len != 0 && (t % print_freq > print_freq-short_len || t % print_freq == 0)){
avg_print_short_e+=energy/short_len;
avg_print_short_a+=alpha/short_len;
avg_print_short_en+=enstrophy/short_len;
}
if(t % print_freq ==0 && t>0){
// compute averages
if (t > running_avg_window) {
avg_e=save_print_short_e[nr_print_in_window]*((double)short_len/running_avg_window);
avg_a=save_print_short_a[nr_print_in_window]*((double)short_len/running_avg_window);
avg_en=save_print_short_en[nr_print_in_window]*((double)short_len/running_avg_window);
for(i=0;i<nr_print_in_window;i++){
// prevent warnings: these will be initialized as long as t > running_avg_window
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
avg_e+=save_print_e[i]*((double)print_freq/running_avg_window);
avg_a+=save_print_a[i]*((double)print_freq/running_avg_window);
avg_en+=save_print_en[i]*((double)print_freq/running_avg_window);
#pragma GCC diagnostic pop
}
}
// memorize averages for future use
// need to keep track of case ==0, otherwise nr_print_in_window-1 is not an unsigned int
if(nr_print_in_window>0){
for(i=0;i<nr_print_in_window-1;i++){
save_print_e[i+1]=save_print_e[i];
save_print_a[i+1]=save_print_a[i];
save_print_en[i+1]=save_print_en[i];
}
}
for(i=0;i<nr_print_in_window;i++){
save_print_short_e[i+1]=save_print_short_e[i];
save_print_short_a[i+1]=save_print_short_a[i];
save_print_short_en[i+1]=save_print_short_en[i];
}
save_print_e[0]=avg_print_e;
save_print_a[0]=avg_print_a;
save_print_en[0]=avg_print_en;
save_print_short_e[0]=avg_print_short_e;
save_print_short_a[0]=avg_print_short_a;
save_print_short_en[0]=avg_print_short_en;
}
}
if(t>0 && t%print_freq==0){
if(t>running_avg_window && t%print_freq==0){
fprintf(stderr,"%d % .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",t,t*delta, avg_a, avg_e, avg_en, alpha, energy, enstrophy);
printf("%8d % .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",t,t*delta, avg_a, avg_e, avg_en, alpha, energy, enstrophy);
}
}
if(running_avg_window!=0){
free(save_print_e);
free(save_print_a);
free(save_print_en);
free(save_print_short_e);
free(save_print_short_a);
free(save_print_short_en);
}
// save final entry to savefile
write_u(u, K1, K2, savefile);