Ian Jauslin

Ian Jauslin

Home Publications Talks Software Teaching Contact
summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/main.c24
-rw-r--r--src/navier-stokes.c121
-rw-r--r--src/navier-stokes.h2
3 files changed, 23 insertions, 124 deletions
diff --git a/src/main.c b/src/main.c
index d7aa8fd..cc1ad32 100644
--- a/src/main.c
+++ b/src/main.c
@@ -28,7 +28,6 @@ typedef struct nstrophy_parameters {
double delta;
double L;
unsigned int print_freq;
- unsigned int avg_window;
int seed;
unsigned int starting_time;
unsigned int driving;
@@ -44,7 +43,7 @@ int print_params(nstrophy_parameters parameters, char* initfile_str, FILE* file)
// read command line arguments
int read_args(int argc, const char* argv[], char** params, unsigned int* command, unsigned int* nthreads, char** savefile_str);
int read_params(char* param_str, nstrophy_parameters* parameters, char** initfile_str);
-int set_parameter(char* lhs, char* rhs, nstrophy_parameters* parameters, bool* setN1, bool* setN2, bool* setavg_window, char** initfile_str);
+int set_parameter(char* lhs, char* rhs, nstrophy_parameters* parameters, bool* setN1, bool* setN2, char** initfile_str);
// set driving force
_Complex double* set_driving(nstrophy_parameters parameters);
@@ -149,7 +148,7 @@ int main (
else if(command==COMMAND_EEA){
// register signal handler to handle aborts
signal(SIGINT, sig_handler);
- 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, (char*)argv[0], param_str, savefile_str);
+ 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, (char*)argv[0], param_str, savefile_str);
}
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);
@@ -323,8 +322,6 @@ 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;
@@ -365,7 +362,7 @@ int read_params(
break;
case ';':
//set parameter
- ret=set_parameter(buffer_lhs, buffer_rhs, parameters, &setN1, &setN2, &setavg_window, initfile_str);
+ ret=set_parameter(buffer_lhs, buffer_rhs, parameters, &setN1, &setN2, initfile_str);
if(ret<0){
return ret;
}
@@ -392,7 +389,7 @@ int read_params(
// set last param
if (*param_str!='\0'){
- ret=set_parameter(buffer_lhs, buffer_rhs, parameters, &setN1, &setN2, &setavg_window, initfile_str);
+ ret=set_parameter(buffer_lhs, buffer_rhs, parameters, &setN1, &setN2, initfile_str);
if(ret<0){
return ret;
}
@@ -410,10 +407,6 @@ 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);
}
@@ -426,7 +419,6 @@ int set_parameter(
nstrophy_parameters* parameters,
bool* setN1,
bool* setN2,
- bool* setavg_window,
char** initfile_str
){
int ret;
@@ -533,14 +525,6 @@ 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 40fb9aa..b9a3979 100644
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -91,7 +91,6 @@ 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,
@@ -106,11 +105,7 @@ 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;
@@ -126,27 +121,8 @@ 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);
- }
+ // special first case when starting_time is not a multiple of print_freq
+ unsigned int first_box = print_freq - (starting_time % print_freq);
// iterate
for(t=starting_time;t<starting_time+nsteps;t++){
@@ -160,78 +136,26 @@ int eea(
enstrophy=compute_enstrophy(u, K1, K2, L);
// running average
- // 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;
- }
+ // reset averages
+ if(t % print_freq == 1){
+ avg_e=0;
+ avg_a=0;
+ avg_en=0;
}
- 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>starting_time){
- // compute averages
- if (t > running_avg_window + starting_time) {
- 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;
- }
+ // compute average
+ // different computationin first box if starting_time is not a multiple of print_freq
+ if(t < starting_time + first_box){
+ avg_e+=energy/first_box;
+ avg_a+=alpha/first_box;
+ avg_en+=enstrophy/first_box;
+ } else {
+ avg_e+=energy/print_freq;
+ avg_a+=alpha/print_freq;
+ avg_en+=enstrophy/print_freq;
}
- if(t>running_avg_window + starting_time && t%print_freq==0){
+ if(t>starting_time && 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);
}
@@ -272,15 +196,6 @@ int eea(
write_u_bin(u, K1, K2, savefile);
}
- 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);
- }
-
ns_free_tmps(u, tmp1, tmp2, tmp3, fft1, fft2, ifft);
return(0);
}
diff --git a/src/navier-stokes.h b/src/navier-stokes.h
index 66178b9..ec304a4 100644
--- a/src/navier-stokes.h
+++ b/src/navier-stokes.h
@@ -20,7 +20,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 running_avg_window, unsigned int starting_time, unsigned int nthreads, FILE* savefile, char* cmd_string, char* params_string, char* savefile_string);
+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, char* cmd_string, char* params_string, char* savefile_string);
// 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);
generated by cgit v1.2.3-70-g09d2 (git 2.51.1) at 2025-10-31 13:19:03 +0000