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

Sort Lyapunov exponents and print more

Browse Source
This commit is contained in:
Ian Jauslin 2024-02-20 09:32:32 -05:00
parent bf8e4b728e
commit c284587105
2 changed files with 30 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
src/lyapunov.c
View File

@ -128,16 +128,29 @@ int lyapunov(
// extract eigenvalues (diagonal elements of Du_prod)
for(i=0; i<MATSIZE; i++){
lyapunov[i]=log(cabs(Du_prod[i*MATSIZE+i]))/(time-prevtime);
// add to average
}
// sort lyapunov exponents
qsort(lyapunov, MATSIZE, sizeof(double), compare_double);
// average lyapunov
for(i=0; i<MATSIZE; i++){
// exclude inf
if(! isinf(lyapunov[i])){
lyapunov_avg[i]=lyapunov_avg[i]*prevtime/time+lyapunov[i]*(time-prevtime)/time;
}
}
// print
fprintf(stderr,"% .15e\n",time);
for(i=0; i<MATSIZE; i++){
printf("% .15e % .15e % .15e\n",time, lyapunov[i], lyapunov_avg[i]);
}
printf("\n");
fprintf(stderr,"% .15e",time);
// print largest and smallest lyapunov exponent
if(MATSIZE>0){
fprintf(stderr," % .15e % .15e\n", lyapunov[0], lyapunov[MATSIZE-1]);
}
// set Du_prod to Q:
LAPACKE_zungrq(LAPACK_COL_MAJOR, MATSIZE, MATSIZE, MATSIZE, Du_prod, MATSIZE, tmp2);
@ -151,6 +164,17 @@ int lyapunov(
return(0);
}
// comparison function for qsort
int compare_double(const void* x, const void* y) {
if (*(const double*)x<*(const double*)y) {
return(-1);
} else if (*(const double*)x>*(const double*)y) {
return(+1);
} else {
return(0);
}
}
// Jacobian of u_n -> u_{n+1}
int ns_D_step(
_Complex double* out,

3
src/lyapunov.h
View File

@ -22,6 +22,9 @@ limitations under the License.
// compute Lyapunov exponents
int lyapunov( int K1, int K2, int N1, int N2, double final_time, double lyapunov_reset, double nu, double D_epsilon, double delta, double L, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_norm, _Complex double* u0, _Complex double* g, bool irreversible, bool keep_en_cst, double target_en, unsigned int algorithm, double starting_time, unsigned int nthreads);
// comparison function for qsort
int compare_double(const void* x, const void* y);
// Jacobian of step
int ns_D_step( _Complex double* out, _Complex double* un, _Complex double* un_next, int K1, int K2, int N1, int N2, double nu, double epsilon, double delta, unsigned int algorithm, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_norm, double L, _Complex double* g, double time, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double* tmp1, _Complex double* tmp2, _Complex double* tmp3, _Complex double* tmp4, _Complex double* tmp5, _Complex double* tmp6, _Complex double* tmp7, _Complex double* tmp8, bool irreversible, bool keep_en_cst, double target_en);