From 59100a471f90b46023741c16fffb69f72e6b8edd Mon Sep 17 00:00:00 2001 From: Ian Jauslin Date: Wed, 12 Apr 2023 14:37:02 -0400 Subject: Running average --- src/main.c | 26 +++++++++--- src/navier-stokes.c | 112 +++++++++++++++++++++++++++++++++++++++++++++++++--- src/navier-stokes.h | 2 +- 3 files changed, 129 insertions(+), 11 deletions(-) (limited to 'src') diff --git a/src/main.c b/src/main.c index 9962c77..cbae30e 100644 --- a/src/main.c +++ b/src/main.c @@ -25,6 +25,7 @@ typedef struct nstrophy_parameters { double delta; double L; unsigned int print_freq; + unsigned int avg_window; int seed; unsigned int starting_time; } nstrophy_parameters; @@ -37,7 +38,7 @@ int print_params(nstrophy_parameters parameters, unsigned int driving, unsigned // read command line arguments int read_args(int argc, const char* argv[], char** params, unsigned int* driving_force, unsigned int* command, unsigned int* init, unsigned int* nthreads, char** savefile_str, char** initfile_str); int read_params(char* param_str, nstrophy_parameters* parameters); -int set_parameter(char* lhs, char* rhs, nstrophy_parameters* parameters, bool* setN1, bool* setN2); +int set_parameter(char* lhs, char* rhs, nstrophy_parameters* parameters, bool* setN1, bool* setN2, bool* setavg_window); // set driving force _Complex double* set_driving(unsigned int driving, nstrophy_parameters parameters); @@ -125,7 +126,7 @@ int main ( uk(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.irreversible, parameters.print_freq, parameters.starting_time, nthreads, savefile); } else if(command==COMMAND_EEA){ - eea(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.irreversible, parameters.print_freq, parameters.starting_time, nthreads, savefile); + eea(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.irreversible, parameters.print_freq, parameters.avg_window, parameters.starting_time, nthreads, savefile); } else if(command==COMMAND_QUIET){ quiet(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.irreversible, nthreads, savefile); @@ -343,6 +344,8 @@ int read_params( // whether N was set explicitly bool setN1=false; bool setN2=false; + // whether avg_window was set explicitly + bool setavg_window=false; // whether lhs (false is rhs) bool lhs=true; @@ -380,7 +383,7 @@ int read_params( break; case ';': //set parameter - ret=set_parameter(buffer_lhs, buffer_rhs, parameters, &setN1, &setN2); + ret=set_parameter(buffer_lhs, buffer_rhs, parameters, &setN1, &setN2, &setavg_window); if(ret<0){ return ret; } @@ -407,7 +410,7 @@ int read_params( // set last param if (*param_str!='\0'){ - ret=set_parameter(buffer_lhs, buffer_rhs, parameters, &setN1, &setN2); + ret=set_parameter(buffer_lhs, buffer_rhs, parameters, &setN1, &setN2, &setavg_window); if(ret<0){ return ret; } @@ -425,6 +428,10 @@ int read_params( if (!setN2){ parameters->N2=smallest_pow2(3*(parameters->K2)); } + // if avg_window is not set explicitly, set it to print_freq + if (!setavg_window){ + parameters->avg_window=parameters->print_freq; + } return(0); } @@ -436,7 +443,8 @@ int set_parameter( char* rhs, nstrophy_parameters* parameters, bool* setN1, - bool* setN2 + bool* setN2, + bool* setavg_window ){ int ret; @@ -542,6 +550,14 @@ int set_parameter( return(-1); } } + else if (strcmp(lhs,"avg_window")==0){ + ret=sscanf(rhs,"%u",&(parameters->avg_window)); + if(ret!=1){ + fprintf(stderr, "error: parameter 'avg_window' should be an integer\n got '%s'\n",rhs); + return(-1); + } + *setavg_window=true; + } else if (strcmp(lhs,"random_seed")==0){ ret=sscanf(rhs,"%d",&(parameters->seed)); if(ret!=1){ diff --git a/src/navier-stokes.c b/src/navier-stokes.c index 19e2f69..5d71d73 100644 --- a/src/navier-stokes.c +++ b/src/navier-stokes.c @@ -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;t0){ - 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 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;i0 && 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); diff --git a/src/navier-stokes.h b/src/navier-stokes.h index d100f44..a8dda69 100644 --- a/src/navier-stokes.h +++ b/src/navier-stokes.h @@ -17,7 +17,7 @@ typedef struct fft_vects { int uk( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int print_freq, unsigned int starting_time, unsigned int nthreadsl, FILE* savefile); // compute energy, enstrophy and alpha -int eea( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int print_freq, unsigned int starting_time, unsigned int nthreads, FILE* savefile); +int eea( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int print_freq, unsigned int running_avg_window, unsigned int starting_time, unsigned int nthreads, FILE* savefile); // 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, double L, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int nthreads, FILE* savefile); -- cgit v1.2.3-54-g00ecf