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

Fix interruption of lyapunov

Browse Source
This commit is contained in:
Ian Jauslin
2025-02-01 11:53:19 -05:00
parent d1992eca70
commit 03c2d1b02a
4 changed files with 40 additions and 19 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/io.c
View File

@@ -58,7 +58,7 @@ int write_vec2_bin(double* vec, int K1, int K2, FILE* file){
return 0; 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; return 0;
} }

31
src/lyapunov.c
View File

@@ -52,6 +52,8 @@ int lyapunov(
unsigned int nthreads, unsigned int nthreads,
double* flow0, double* flow0,
double* lyapunov_avg0, 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* savefile,
FILE* utfile, FILE* utfile,
// for interrupt recovery // for interrupt recovery
@@ -105,9 +107,6 @@ int lyapunov(
double step=delta; double step=delta;
double next_step=step; double next_step=step;
// save times at which Lyapunov exponents are computed
double prevtime=starting_time;
// iterate // iterate
time=starting_time; time=starting_time;
while(final_time<0 || time<final_time){ while(final_time<0 || time<final_time){
@@ -155,8 +154,8 @@ int lyapunov(
// average lyapunov // average lyapunov
for(i=0; i<MATSIZE; i++){ for(i=0; i<MATSIZE; i++){
// exclude inf // exclude inf
if((! isinf(lyapunov[i])) && (time>starting_time)){ if((! isinf(lyapunov[i])) && (time>lyapunov_startingtime)){
lyapunov_avg[i]=lyapunov_avg[i]*(prevtime-starting_time)/(time-starting_time)+lyapunov[i]*(time-prevtime)/(time-starting_time); 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 // reset prevtime
prevtime=time; prevtime=time;
}
// catch abort signal // catch abort signal
if (g_abort){ if (g_abort){
// print u to stderr if no savefile // print u to stderr if no savefile
if (savefile==NULL){ if (savefile==NULL){
savefile=stderr; savefile=stderr;
}
break;
} }
break;
} }
} }
if(savefile!=NULL){ 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 // save final u to utfile in txt format
@@ -781,6 +780,8 @@ int lyapunov_save_state(
double* flow, double* flow,
_Complex double* u, _Complex double* u,
double* lyapunov_avg, double* lyapunov_avg,
double prevtime,
double lyapunov_startingtime,
FILE* savefile, FILE* savefile,
int K1, int K1,
int K2, int K2,
@@ -801,6 +802,10 @@ int lyapunov_save_state(
if(savefile!=stderr && savefile!=stdout){ if(savefile!=stderr && savefile!=stdout){
// save tangent flow // save tangent flow
write_mat2_bin(flow,K1,K2,savefile); 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 // save lyapunov averages
write_vec2_bin(lyapunov_avg,K1,K2,savefile); write_vec2_bin(lyapunov_avg,K1,K2,savefile);
} }

4
src/lyapunov.h
View File

@@ -20,7 +20,7 @@ limitations under the License.
#include "navier-stokes.h" #include "navier-stokes.h"
// compute Lyapunov exponents // 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 // comparison function for qsort
int compare_double(const void* x, const void* y); 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); 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 // 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 #endif

22
src/main.c
View File

@@ -85,7 +85,7 @@ _Complex double* set_driving(nstrophy_parameters parameters);
// set initial condition // set initial condition
_Complex double* set_init(nstrophy_parameters* parameters); _Complex double* set_init(nstrophy_parameters* parameters);
// set initial tangent flow for lyapunov exponents // 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 // signal handler
void sig_handler (int signo); void sig_handler (int signo);
@@ -121,6 +121,8 @@ int main (
_Complex double *g; _Complex double *g;
double* lyapunov_flow0; double* lyapunov_flow0;
double* lyapunov_avg0; double* lyapunov_avg0;
double lyapunov_prevtime;
double lyapunov_startingtime;
dstring savefile_str; dstring savefile_str;
dstring utfile_str; dstring utfile_str;
dstring initfile_str; dstring initfile_str;
@@ -232,7 +234,7 @@ int main (
u0=set_init(&parameters); u0=set_init(&parameters);
// set initial condition for the lyapunov flow // set initial condition for the lyapunov flow
if (command==COMMAND_LYAPUNOV){ 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 // 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); 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){ 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(&param_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(&param_str), dstring_to_str_noinit(&savefile_str), dstring_to_str_noinit(&utfile_str));
} }
else if(command==0){ else if(command==0){
fprintf(stderr, "error: no command specified\n"); fprintf(stderr, "error: no command specified\n");
@@ -1170,6 +1175,8 @@ _Complex double* set_init(
int set_lyapunov_flow_init( int set_lyapunov_flow_init(
double** lyapunov_flow0, double** lyapunov_flow0,
double** lyapunov_avg0, double** lyapunov_avg0,
double* lyapunov_prevtime,
double* lyapunov_startingtime,
nstrophy_parameters parameters 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_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){ if(parameters.init_flow_file){
// read flow // read flow
read_mat2_bin(*lyapunov_flow0, parameters.K1, parameters.K2, parameters.initfile); 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 averages
read_vec2_bin(*lyapunov_avg0, parameters.K1, parameters.K2, parameters.initfile); read_vec2_bin(*lyapunov_avg0, parameters.K1, parameters.K2, parameters.initfile);
} else { } 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++){ for(i=0;i<2*(parameters.K1*(2*parameters.K2+1)+parameters.K2);i++){
(*lyapunov_avg0)[i]=0; (*lyapunov_avg0)[i]=0;
} }