diff --git a/src/io.c b/src/io.c index fb05327..6932a8f 100644 --- a/src/io.c +++ b/src/io.c @@ -58,7 +58,7 @@ int write_vec2_bin(double* vec, int K1, int K2, FILE* file){ return 0; } - fwrite(vec, sizeof(double), 2*K1*(2*K2+1)+K2, file); + fwrite(vec, sizeof(double), 2*(K1*(2*K2+1)+K2), file); return 0; } diff --git a/src/lyapunov.c b/src/lyapunov.c index 99c6ffc..bab013d 100644 --- a/src/lyapunov.c +++ b/src/lyapunov.c @@ -52,6 +52,8 @@ int lyapunov( unsigned int nthreads, double* flow0, double* lyapunov_avg0, + double prevtime, // the previous time at which a QR decomposition was performed + double lyapunov_startingtime, // the time at which the lyapunov exponent computation was started (useful in interrupted computation) FILE* savefile, FILE* utfile, // for interrupt recovery @@ -105,9 +107,6 @@ int lyapunov( double step=delta; double next_step=step; - // save times at which Lyapunov exponents are computed - double prevtime=starting_time; - // iterate time=starting_time; while(final_time<0 || timestarting_time)){ - lyapunov_avg[i]=lyapunov_avg[i]*(prevtime-starting_time)/(time-starting_time)+lyapunov[i]*(time-prevtime)/(time-starting_time); + if((! isinf(lyapunov[i])) && (time>lyapunov_startingtime)){ + lyapunov_avg[i]=lyapunov_avg[i]*(prevtime-lyapunov_startingtime)/(time-lyapunov_startingtime)+lyapunov[i]*(time-prevtime)/(time-lyapunov_startingtime); } } @@ -176,20 +175,20 @@ int lyapunov( // reset prevtime prevtime=time; + } - // catch abort signal - if (g_abort){ - // print u to stderr if no savefile - if (savefile==NULL){ - savefile=stderr; - } - break; + // catch abort signal + if (g_abort){ + // print u to stderr if no savefile + if (savefile==NULL){ + savefile=stderr; } + break; } } if(savefile!=NULL){ - lyapunov_save_state(flow, u, lyapunov_avg, savefile, K1, K2, cmd_string, params_string, savefile_string, utfile_string, utfile, COMMAND_LYAPUNOV, algorithm, step, time, nthreads); + lyapunov_save_state(flow, u, lyapunov_avg, prevtime, lyapunov_startingtime, savefile, K1, K2, cmd_string, params_string, savefile_string, utfile_string, utfile, COMMAND_LYAPUNOV, algorithm, step, time, nthreads); } // save final u to utfile in txt format @@ -781,6 +780,8 @@ int lyapunov_save_state( double* flow, _Complex double* u, double* lyapunov_avg, + double prevtime, + double lyapunov_startingtime, FILE* savefile, int K1, int K2, @@ -801,6 +802,10 @@ int lyapunov_save_state( if(savefile!=stderr && savefile!=stdout){ // save tangent flow write_mat2_bin(flow,K1,K2,savefile); + // save time of QR decomposition + fwrite(&prevtime, sizeof(double), 1, savefile); + // save time at which the lyapunov computation started + fwrite(&lyapunov_startingtime, sizeof(double), 1, savefile); // save lyapunov averages write_vec2_bin(lyapunov_avg,K1,K2,savefile); } diff --git a/src/lyapunov.h b/src/lyapunov.h index 38b0714..c37ebbf 100644 --- a/src/lyapunov.h +++ b/src/lyapunov.h @@ -20,7 +20,7 @@ limitations under the License. #include "navier-stokes.h" // compute Lyapunov exponents -int lyapunov( int K1, int K2, int N1, int N2, double final_time, double lyapunov_reset, unsigned int lyapunov_trigger, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_cost, _Complex double* u0, _Complex double* g, bool irreversible, bool keep_en_cst, double target_en, unsigned int algorithm, unsigned int algorithm_lyapunov, double starting_time, unsigned int nthreads, double* flow0, double* lyapunov_avg0, FILE* savefile, FILE* utfile, const char* cmd_string, const char* params_string, const char* savefile_string, const char* utfile_string); +int lyapunov( int K1, int K2, int N1, int N2, double final_time, double lyapunov_reset, unsigned int lyapunov_trigger, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_cost, _Complex double* u0, _Complex double* g, bool irreversible, bool keep_en_cst, double target_en, unsigned int algorithm, unsigned int algorithm_lyapunov, double starting_time, unsigned int nthreads, double* flow0, double* lyapunov_avg0, double prevtime, double lyapunov_startingtime, FILE* savefile, FILE* utfile, const char* cmd_string, const char* params_string, const char* savefile_string, const char* utfile_string); // comparison function for qsort int compare_double(const void* x, const void* y); @@ -57,7 +57,7 @@ int lyapunov_init_tmps( double ** lyapunov, double ** lyapunov_avg, double ** fl int lyapunov_free_tmps( double * lyapunov, double * lyapunov_avg, double * flow, _Complex double * u, double * tmp1, double * tmp2, double * tmp3, _Complex double * tmp4, _Complex double * tmp5, _Complex double * tmp6, _Complex double * tmp7, _Complex double * tmp8, _Complex double * tmp9, _Complex double * tmp10, double * tmp11, fft_vect fftu1, fft_vect fftu2, fft_vect fftu3, fft_vect fftu4, fft_vect fft1, fft_vect ifft, unsigned int algorithm, unsigned int algorithm_lyapunov); // save state of lyapunov computation -int lyapunov_save_state( double* flow, _Complex double* u, double* lyapunov_avg, FILE* savefile, int K1, int K2, const char* cmd_string, const char* params_string, const char* savefile_string, const char* utfile_string, FILE* utfile, unsigned int command, unsigned int algorithm, double step, double time, unsigned int nthreads); +int lyapunov_save_state( double* flow, _Complex double* u, double* lyapunov_avg, double prevtime, double lyapunov_startingtime, FILE* savefile, int K1, int K2, const char* cmd_string, const char* params_string, const char* savefile_string, const char* utfile_string, FILE* utfile, unsigned int command, unsigned int algorithm, double step, double time, unsigned int nthreads); #endif diff --git a/src/main.c b/src/main.c index 74c6f7d..2d5651a 100644 --- a/src/main.c +++ b/src/main.c @@ -85,7 +85,7 @@ _Complex double* set_driving(nstrophy_parameters parameters); // set initial condition _Complex double* set_init(nstrophy_parameters* parameters); // set initial tangent flow for lyapunov exponents -int set_lyapunov_flow_init( double** lyapunov_flow0, double** lyapunov_avg0, nstrophy_parameters parameters); +int set_lyapunov_flow_init( double** lyapunov_flow0, double** lyapunov_avg0, double* lyapunov_prevtime, double* lyapunov_startingtime, nstrophy_parameters parameters); // signal handler void sig_handler (int signo); @@ -121,6 +121,8 @@ int main ( _Complex double *g; double* lyapunov_flow0; double* lyapunov_avg0; + double lyapunov_prevtime; + double lyapunov_startingtime; dstring savefile_str; dstring utfile_str; dstring initfile_str; @@ -232,7 +234,7 @@ int main ( u0=set_init(¶meters); // set initial condition for the lyapunov flow if (command==COMMAND_LYAPUNOV){ - set_lyapunov_flow_init(&lyapunov_flow0, &lyapunov_avg0, parameters); + set_lyapunov_flow_init(&lyapunov_flow0, &lyapunov_avg0, &lyapunov_prevtime, &lyapunov_startingtime, parameters); } // if init_enstrophy is not set in the parameters, then compute it from the initial condition @@ -315,7 +317,10 @@ int main ( quiet(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_cost, parameters.starting_time, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_enstrophy, parameters.algorithm, nthreads, savefile); } else if(command==COMMAND_LYAPUNOV){ - lyapunov(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.lyapunov_reset, parameters.lyapunov_trigger, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_cost, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_enstrophy, parameters.algorithm, parameters.algorithm_lyapunov, parameters.starting_time, nthreads, lyapunov_flow0, lyapunov_avg0, savefile, utfile, (char*)argv[0], dstring_to_str_noinit(¶m_str), dstring_to_str_noinit(&savefile_str), dstring_to_str_noinit(&utfile_str)); + // register signal handler to handle aborts + signal(SIGINT, sig_handler); + signal(SIGTERM, sig_handler); + lyapunov(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.lyapunov_reset, parameters.lyapunov_trigger, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_cost, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_enstrophy, parameters.algorithm, parameters.algorithm_lyapunov, parameters.starting_time, nthreads, lyapunov_flow0, lyapunov_avg0, lyapunov_prevtime, lyapunov_startingtime, savefile, utfile, (char*)argv[0], dstring_to_str_noinit(¶m_str), dstring_to_str_noinit(&savefile_str), dstring_to_str_noinit(&utfile_str)); } else if(command==0){ fprintf(stderr, "error: no command specified\n"); @@ -1170,6 +1175,8 @@ _Complex double* set_init( int set_lyapunov_flow_init( double** lyapunov_flow0, double** lyapunov_avg0, + double* lyapunov_prevtime, + double* lyapunov_startingtime, nstrophy_parameters parameters ){ *lyapunov_flow0=calloc(4*(parameters.K1*(2*parameters.K2+1)+parameters.K2)*(parameters.K1*(2*parameters.K2+1)+parameters.K2),sizeof(double)); @@ -1178,6 +1185,10 @@ int set_lyapunov_flow_init( if(parameters.init_flow_file){ // read flow read_mat2_bin(*lyapunov_flow0, parameters.K1, parameters.K2, parameters.initfile); + // read time of last QR decomposition + fread(lyapunov_prevtime, sizeof(double), 1, parameters.initfile); + // read time at which lyapunov computation was started + fread(lyapunov_startingtime, sizeof(double), 1, parameters.initfile); // read averages read_vec2_bin(*lyapunov_avg0, parameters.K1, parameters.K2, parameters.initfile); } else { @@ -1192,6 +1203,11 @@ int set_lyapunov_flow_init( } } } + // init prevtime + *lyapunov_prevtime=parameters.starting_time; + // init starting_time + *lyapunov_startingtime=parameters.starting_time; + // init averages for(i=0;i<2*(parameters.K1*(2*parameters.K2+1)+parameters.K2);i++){ (*lyapunov_avg0)[i]=0; }