From 06e5b0e0da9754468ebdaa72f9534b6a6669ccab Mon Sep 17 00:00:00 2001 From: Ian Jauslin Date: Fri, 31 Jan 2025 21:58:54 -0500 Subject: [PATCH] Allow the Lyapunov computation to be interrupted --- src/io.c | 74 ++++++++++++++++++++++++++++++++++++++++++++++---- src/io.h | 11 ++++++++ src/lyapunov.c | 72 +++++++++++++++++++++++++++++++++++++++++++----- src/lyapunov.h | 5 +++- src/main.c | 70 +++++++++++++++++++++++++++++++++++++++++++++-- 5 files changed, 217 insertions(+), 15 deletions(-) diff --git a/src/io.c b/src/io.c index 127e1bc..fb05327 100644 --- a/src/io.c +++ b/src/io.c @@ -51,6 +51,30 @@ int write_vec_bin(_Complex double* vec, int K1, int K2, FILE* file){ return 0; } +// write complex vector (stored as 2 doubles) indexed by k1,k2 to file in binary format +int write_vec2_bin(double* vec, int K1, int K2, FILE* file){ + // do nothing if there is no file + if(file==NULL){ + return 0; + } + + fwrite(vec, sizeof(double), 2*K1*(2*K2+1)+K2, file); + + return 0; +} + +// write complex matrix (stored as 2 doubles) indexed by k1,k2 to file in binary format +int write_mat2_bin(double* mat, int K1, int K2, FILE* file){ + // do nothing if there is no file + if(file==NULL){ + return 0; + } + + fwrite(mat, sizeof(double), 4*(K1*(2*K2+1)+K2)*(K1*(2*K2+1)+K2), file); + + return 0; +} + // read complex vector indexed by k1,k2 from file int read_vec(_Complex double* out, int K1, int K2, FILE* file){ int kx,ky; @@ -149,15 +173,53 @@ int read_vec(_Complex double* out, int K1, int K2, FILE* file){ // read complex vector indexed by k1,k2 from file in binary format int read_vec_bin(_Complex double* out, int K1, int K2, FILE* file){ - char c; - int ret; - // do nothing if there is no file if(file==NULL){ return 0; } // seek past initial comments + seek_past_initial_comments(file); + + fread(out, sizeof(_Complex double), K1*(2*K2+1)+K2, file); + + return 0; +} + +// read complex vector (represented as 2 doubles) indexed by k1,k2 from file in binary format +int read_vec2_bin(double* out, int K1, int K2, FILE* file){ + if(file==NULL){ + return 0; + } + + // seek past initial comments + seek_past_initial_comments(file); + + fread(out, sizeof(double), 2*(K1*(2*K2+1)+K2), file); + return 0; +} + +// read complex matrix (represented as 2 doubles) indexed by k1,k2 from file in binary format +int read_mat2_bin(double* out, int K1, int K2, FILE* file){ + if(file==NULL){ + return 0; + } + + // seek past initial comments + seek_past_initial_comments(file); + + fread(out, sizeof(double), 4*(K1*(2*K2+1)+K2)*(K1*(2*K2+1)+K2), file); + return 0; +} + +// ignore comments at beginning of file +int seek_past_initial_comments(FILE* file){ + char c; + int ret; + + if(file==NULL){ + return 0; + } while(true){ ret=fscanf(file, "%c", &c); if (ret==1 && c=='#'){ @@ -184,8 +246,6 @@ int read_vec_bin(_Complex double* out, int K1, int K2, FILE* file){ } } - fread(out, sizeof(_Complex double), K1*(2*K2+1)+K2, file); - return 0; } @@ -275,6 +335,10 @@ int save_state( if(savefile==stderr || savefile==stdout){ fprintf(savefile,";starting_time=%.15e", time); } + // instruction to read init flow from file if computation is lyapunov + if(command==COMMAND_LYAPUNOV){ + fprintf(savefile,";init_flow=file"); + } fprintf(savefile,"\""); } diff --git a/src/io.h b/src/io.h index ba17413..e21d56f 100644 --- a/src/io.h +++ b/src/io.h @@ -22,10 +22,21 @@ limitations under the License. // write complex vector indexed by k1,k2 to file int write_vec(_Complex double* u, int K1, int K2, FILE* file); int write_vec_bin(_Complex double* u, int K1, int K2, FILE* file); +// write complex vector (stored as 2 doubles) indexed by k1,k2 to file in binary format +int write_vec2_bin(double* vec, int K1, int K2, FILE* file); +// write complex matrix (stored as 2 doubles) indexed by k1,k2 to file in binary format +int write_mat2_bin(double* mat, int K1, int K2, FILE* file); // read complex vector indexed by k1,k2 from file int read_vec(_Complex double* u, int K1, int K2, FILE* file); int read_vec_bin(_Complex double* u, int K1, int K2, FILE* file); +// read complex vector (represented as 2 doubles) indexed by k1,k2 from file in binary format +int read_vec2_bin(double* out, int K1, int K2, FILE* file); +// read complex matrix (represented as 2 doubles) indexed by k1,k2 from file in binary format +int read_mat2_bin(double* out, int K1, int K2, FILE* file); + +// ignore comments at beginning of file +int seek_past_initial_comments(FILE* file); // remove an entry from params string (inplace) int remove_entry(char* param_str, char* entry); diff --git a/src/lyapunov.c b/src/lyapunov.c index d2c48d3..99c6ffc 100644 --- a/src/lyapunov.c +++ b/src/lyapunov.c @@ -16,6 +16,7 @@ limitations under the License. #include "constants.cpp" #include "lyapunov.h" +#include "io.h" #include #include #include @@ -48,7 +49,16 @@ int lyapunov( unsigned int algorithm, unsigned int algorithm_lyapunov, double starting_time, - unsigned int nthreads + unsigned int nthreads, + double* flow0, + double* lyapunov_avg0, + FILE* savefile, + FILE* utfile, + // for interrupt recovery + const char* cmd_string, + const char* params_string, + const char* savefile_string, + const char* utfile_string ){ double* lyapunov; double* lyapunov_avg; @@ -82,15 +92,13 @@ int lyapunov( // copy initial condition copy_u(u, u0, K1, K2); - // initialize flow with identity matrix + + // initialize flow and averages for (i=0;ilength==0){ + if(*command==COMMAND_RESUME && (resumefile_str==NULL || resumefile_str->length==0)){ fprintf(stderr, "error: 'resume' command used, but no resume file\n"); print_usage(); return(-1); @@ -580,6 +602,7 @@ int set_default_params( parameters->algorithm_lyapunov=ALGORITHM_RK4; // default lyapunov_reset will be print_time (set later) for now set target to 0 to indicate it hasn't been set explicitly parameters->lyapunov_trigger=0; + parameters->init_flow_file=false; parameters->print_alpha=false; @@ -977,6 +1000,16 @@ int set_parameter( return(-1); } } + else if (strcmp(lhs,"init_flow")==0){ + if(strcmp(rhs,"file")==0){ + parameters->init_flow_file=true; + } else if (strcmp(rhs,"identity")==0){ + parameters->init_flow_file=false; + } else { + fprintf(stderr, "error: parameter 'init_flow' should be 'file' or 'identity'\n got '%s'\n",rhs); + return(-1); + } + } else{ fprintf(stderr, "error: unrecognized parameter '%s'\n",lhs); return(-1); @@ -1132,3 +1165,36 @@ _Complex double* set_init( return u0; } + +// set initial tangent flow for lyapunov exponents +int set_lyapunov_flow_init( + double** lyapunov_flow0, + double** lyapunov_avg0, + 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)); + *lyapunov_avg0=calloc(2*(parameters.K1*(2*parameters.K2+1)+parameters.K2),sizeof(double)); + + if(parameters.init_flow_file){ + // read flow + read_mat2_bin(*lyapunov_flow0, parameters.K1, parameters.K2, parameters.initfile); + // read averages + read_vec2_bin(*lyapunov_avg0, parameters.K1, parameters.K2, parameters.initfile); + } else { + // init with identity + int i,j; + for (i=0;i<2*(parameters.K1*(2*parameters.K2+1)+parameters.K2);i++){ + for (j=0;j<2*(parameters.K1*(2*parameters.K2+1)+parameters.K2);j++){ + if(i!=j){ + (*lyapunov_flow0)[i*2*(parameters.K1*(2*parameters.K2+1)+parameters.K2)+j]=0.; + } else { + (*lyapunov_flow0)[i*2*(parameters.K1*(2*parameters.K2+1)+parameters.K2)+j]=1.; + } + } + } + for(i=0;i<2*(parameters.K1*(2*parameters.K2+1)+parameters.K2);i++){ + (*lyapunov_avg0)[i]=0; + } + } + return 0; +}