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

Allow the Lyapunov computation to be interrupted

Browse Source
This commit is contained in:
Ian Jauslin
2025-01-31 21:58:54 -05:00
parent 53a0a0ae4c
commit 06e5b0e0da
5 changed files with 217 additions and 15 deletions
Download Patch File Download Diff File Expand all files Collapse all files

70
src/main.c
View File

@@ -28,6 +28,7 @@ limitations under the License.
#include "dstring.h"
#include "init.h"
#include "int_tools.h"
#include "io.h"
#include "lyapunov.h"
#include "navier-stokes.h"
@@ -62,6 +63,7 @@ typedef struct nstrophy_parameters {
FILE* drivingfile;
double lyapunov_reset;
unsigned int lyapunov_trigger;
bool init_flow_file;
bool print_alpha;
} nstrophy_parameters;
@@ -82,6 +84,8 @@ int args_from_file(dstring* params, unsigned int* command, unsigned int* nthread
_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);
// signal handler
void sig_handler (int signo);
@@ -115,6 +119,8 @@ int main (
unsigned int nthreads=1;
_Complex double* u0;
_Complex double *g;
double* lyapunov_flow0;
double* lyapunov_avg0;
dstring savefile_str;
dstring utfile_str;
dstring initfile_str;
@@ -224,6 +230,10 @@ int main (
g=set_driving(parameters);
// set initial condition
u0=set_init(&parameters);
// set initial condition for the lyapunov flow
if (command==COMMAND_LYAPUNOV){
set_lyapunov_flow_init(&lyapunov_flow0, &lyapunov_avg0, parameters);
}
// if init_enstrophy is not set in the parameters, then compute it from the initial condition
if (parameters.init_enstrophy_or_energy!=FIX_ENSTROPHY){
@@ -255,6 +265,10 @@ int main (
dstring_free(drivingfile_str);
free(g);
free(u0);
if (command==COMMAND_LYAPUNOV){
free(lyapunov_flow0);
free(lyapunov_avg0);
}
return(-1);
}
}
@@ -271,6 +285,10 @@ int main (
dstring_free(drivingfile_str);
free(g);
free(u0);
if (command==COMMAND_LYAPUNOV){
free(lyapunov_flow0);
free(lyapunov_avg0);
}
return(-1);
}
}
@@ -297,7 +315,7 @@ 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(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));
}
else if(command==0){
fprintf(stderr, "error: no command specified\n");
@@ -306,6 +324,10 @@ int main (
free(g);
free(u0);
if (command==COMMAND_LYAPUNOV){
free(lyapunov_flow0);
free(lyapunov_avg0);
}
// free strings
dstring_free(param_str);
@@ -540,7 +562,7 @@ int read_args(
}
// check that if the command is 'resume', then resumefile has been set
if(*command==COMMAND_RESUME && resumefile_str->length==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;
}