/* Copyright 2017-2025 Ian Jauslin Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0 Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License. */ #include "constants.cpp" #include "lyapunov.h" #include "io.h" #include #include #include #include #define MATSIZE (2*(K1*(2*K2+1)+K2)) #define USIZE (K1*(2*K2+1)+K2) // 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, 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 const char* cmd_string, const char* params_string, const char* savefile_string, const char* utfile_string ){ double* lyapunov; double* lyapunov_avg; double* flow; _Complex double* u; double time; fft_vect fftu1; fft_vect fftu2; fft_vect fftu3; fft_vect fftu4; fft_vect fft1; fft_vect ifft; 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; int i,j; // period (only useful with LYAPUNOV_TRIGGER_TIME, but compute it anyways: it won't take much time...) // add 0.1 to ensure proper rounding uint64_t n=(uint64_t)((starting_time-fmod(starting_time, lyapunov_reset))/lyapunov_reset+0.1); lyapunov_init_tmps(&lyapunov, &lyapunov_avg, &flow, &u, &tmp1, &tmp2, &tmp3, &tmp4, &tmp5, &tmp6, &tmp7, &tmp8, &tmp9, &tmp10, &tmp11, &fftu1, &fftu2, &fftu3, &fftu4, &fft1, &ifft, K1, K2, N1, N2, nthreads, algorithm, algorithm_lyapunov); // copy initial condition copy_u(u, u0, K1, K2); // initialize flow and averages for (i=0;ilyapunov_reset){ break; } } if(sqrt(norm)>lyapunov_reset){ break; } } } // QR decomposition // Do it also if it is the last step if((lyapunov_trigger==LYAPUNOV_TRIGGER_TIME && time>(n+1)*lyapunov_reset) || (lyapunov_trigger==LYAPUNOV_TRIGGER_SIZE && norm>lyapunov_reset) || time>=final_time){ n++; // QR decomposition // do not touch tmp1, it might be used in the next step LAPACKE_dgeqrf(LAPACK_COL_MAJOR, MATSIZE, MATSIZE, flow, MATSIZE, tmp11); // extract diagonal elements of R (represented as diagonal elements of flow for(i=0; ilyapunov_startingtime)){ lyapunov_avg[i]=lyapunov_avg[i]*(prevtime-lyapunov_startingtime)/(time-lyapunov_startingtime)+lyapunov[i]*(time-prevtime)/(time-lyapunov_startingtime); } } // print for(i=0; i0){ // fprintf(stderr," % .15e % .15e\n", lyapunov[0], lyapunov[MATSIZE-1]); //} // set flow to Q: LAPACKE_dorgqr(LAPACK_COL_MAJOR, MATSIZE, MATSIZE, MATSIZE, flow, MATSIZE, tmp11); // reset prevtime prevtime=time; } // 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, 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 if(utfile!=NULL){ write_vec(u, K1, K2, utfile); } lyapunov_free_tmps(lyapunov, lyapunov_avg, flow, u, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9, tmp10, tmp11, fftu1, fftu2, fftu3, fftu4, fft1, ifft, algorithm, algorithm_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); } } // compute tangent flow int lyapunov_D_step( double* flow, _Complex double* u, int K1, int K2, int N1, int N2, double nu, double delta, unsigned int algorithm, double L, _Complex double* g, fft_vect fftu1, fft_vect fftu2, fft_vect fftu3, fft_vect fftu4, fft_vect fft1, fft_vect ifft, double* tmp1, double* tmp2, double* tmp3, _Complex double* tmp4, bool irreversible ){ int lx,ly; double alpha; int n; // compute fft of u for future use lyapunov_fft_u_T(u,K1,K2,N1,N2,fftu1,fftu2,fftu3,fftu4); // compute T lyapunov_T(N1,N2,fftu1,fftu2,fftu3,fftu4,ifft); // save to vect for(lx=0;lx<=K1;lx++){ for(ly=(lx>0 ? -K2 : 1);ly<=K2;ly++){ tmp4[klookup_sym(lx,ly,K2)]=ifft.fft[klookup(lx,ly,N1,N2)]; } } //compute alpha if (irreversible) { alpha=nu; } else { alpha=compute_alpha(u,K1,K2,g,L); } // loop over columns for(lx=0;lx<=K1;lx++){ for(ly=(lx>0 ? -K2 : 1);ly<=K2;ly++){ for(n=0;n<=1;n++){ // do not use adaptive step algorithms for the tangent flow: it must be locked to the same times as u if(algorithm==ALGORITHM_RK2){ lyapunov_D_step_rk2(flow+(2*klookup_sym(lx,ly,K2)+n)*MATSIZE, u, K1, K2, N1, N2, alpha, delta, L, g, tmp4, fftu1, fftu2, fftu3, fftu4, fft1, ifft, tmp1, tmp2, irreversible); } else if (algorithm==ALGORITHM_RK4) { lyapunov_D_step_rk4(flow+(2*klookup_sym(lx,ly,K2)+n)*MATSIZE, u, K1, K2, N1, N2, alpha, delta, L, g, tmp4, fftu1, fftu2, fftu3, fftu4, fft1, ifft, tmp1, tmp2, tmp3, irreversible); } else { fprintf(stderr,"bug: unknown algorithm for tangent flow: %u, contact ian.jauslin@rutgers.edu\n",algorithm); } } } } return(0); } // RK 4 algorithm int lyapunov_D_step_rk4( double* flow, _Complex double* u, int K1, int K2, int N1, int N2, double nu, double delta, double L, _Complex double* g, _Complex double* T, fft_vect fftu1, fft_vect fftu2, fft_vect fftu3, fft_vect fftu4, fft_vect fft1, fft_vect ifft, double* tmp1, double* tmp2, double* tmp3, bool irreversible ){ int i; // k1 lyapunov_D_rhs(tmp1, flow, u, K1, K2, N1, N2, nu, L, g, T, fftu1, fftu2, fftu3, fftu4, fft1, ifft, irreversible); // add to output for(i=0;i0 ? -K2 : 1);ky<=K2;ky++){ // real part out[2*klookup_sym(kx,ky,K2)]=-4*M_PI*M_PI/L/L*alpha*(kx*kx+ky*ky)*flow[2*klookup_sym(kx,ky,K2)]+4*M_PI*M_PI/L/L/sqrt(kx*kx+ky*ky)*__real__ ifft.fft[klookup(kx,ky,N1,N2)]; // imaginary part out[2*klookup_sym(kx,ky,K2)+1]=-4*M_PI*M_PI/L/L*alpha*(kx*kx+ky*ky)*flow[2*klookup_sym(kx,ky,K2)+1]+4*M_PI*M_PI/L/L/sqrt(kx*kx+ky*ky)*__imag__ ifft.fft[klookup(kx,ky,N1,N2)]; } } if(!irreversible){ double Dalpha=lyapunov_D_alpha(flow,u,K1,K2,N1,N2,alpha,L,g,T,ifft.fft); for(kx=0;kx<=K1;kx++){ for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){ // real part out[2*klookup_sym(kx,ky,K2)]+=4*M_PI*M_PI/L/L*(kx*kx+ky*ky)*Dalpha*__real__ u[klookup_sym(kx,ky,K2)]; // imaginary part out[2*klookup_sym(kx,ky,K2)+1]+=4*M_PI*M_PI/L/L*(kx*kx+ky*ky)*Dalpha*__imag__ u[klookup_sym(kx,ky,K2)]; } } } return(0); } // Compute T from the already computed fourier transforms of u int lyapunov_T( int N1, int N2, fft_vect fftu1, fft_vect fftu2, fft_vect fftu3, fft_vect fftu4, fft_vect ifft ){ int i; for(i=0;i0 ? -K2 : 1);ky<=K2;ky++){ num+=L*L/4/M_PI/M_PI*(kx*kx+ky*ky)*(flow[2*klookup_sym(kx,ky,K2)]*(__real__ g[klookup_sym(kx,ky,K2)])+flow[2*klookup_sym(kx,ky,K2)+1]*(__imag__ g[klookup_sym(kx,ky,K2)]))// +sqrt(kx*kx+ky*ky)*(flow[2*klookup_sym(kx,ky,K2)]*(__real__ T[klookup_sym(kx,ky,K2)])+flow[2*klookup_sym(kx,ky,K2)+1]*(__imag__ T[klookup_sym(kx,ky,K2)])// // recall that DT is ifft.fft, so is of size N1xN2 +__real__(conj(u[klookup_sym(kx,ky,K2)])*DT[klookup(kx,ky,N1,N2)]))// -2*alpha*(kx*kx+ky*ky)*(kx*kx+ky*ky)*(flow[2*klookup_sym(kx,ky,K2)]*(__real__ u[klookup_sym(kx,ky,K2)])+flow[2*klookup_sym(kx,ky,K2)+1]*(__imag__ u[klookup_sym(kx,ky,K2)])); denom+=(kx*kx+ky*ky)*(kx*kx+ky*ky)*__real__(u[klookup_sym(kx,ky,K2)]*conj(u[klookup_sym(kx,ky,K2)])); } } return num/denom; } // get delta_{kx,ky} from a vector delta in which only the values for kx>=0 are stored, as separate real part and imaginary part _Complex double lyapunov_delta_getval_sym( double* delta, int kx, int ky, int K2 ){ if(kx>0 || (kx==0 && ky>0)){ return delta[2*klookup_sym(kx,ky,K2)]+delta[2*klookup_sym(kx,ky,K2)+1]*_Complex_I; } else if(kx==0 && ky==0){ return 0; } else { return delta[2*klookup_sym(-kx,-ky,K2)]-delta[2*klookup_sym(-kx,-ky,K2)+1]*_Complex_I; } } // compute ffts of u for future use in DT int lyapunov_fft_u_T( _Complex double* u, int K1, int K2, int N1, int N2, fft_vect fftu1, fft_vect fftu2, fft_vect fftu3, fft_vect fftu4 ){ int kx,ky; int i; // F(px/|p|*u)*F(qy*|q|*u) // init to 0 for(i=0; ifft=fftw_malloc(sizeof(fftw_complex)*N1*N2); fftu2->fft_plan=fftw_plan_dft_2d(N1,N2, fftu2->fft, fftu2->fft, FFTW_FORWARD, FFTW_MEASURE); fftu3->fft=fftw_malloc(sizeof(fftw_complex)*N1*N2); fftu3->fft_plan=fftw_plan_dft_2d(N1,N2, fftu3->fft, fftu3->fft, FFTW_FORWARD, FFTW_MEASURE); fftu4->fft=fftw_malloc(sizeof(fftw_complex)*N1*N2); fftu4->fft_plan=fftw_plan_dft_2d(N1,N2, fftu4->fft, fftu4->fft, FFTW_FORWARD, FFTW_MEASURE); *lyapunov=calloc(MATSIZE,sizeof(double)); *lyapunov_avg=calloc(MATSIZE,sizeof(double)); *flow=calloc(MATSIZE*MATSIZE,sizeof(double)); return(0); } // release vectors 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 ){ free(flow); free(lyapunov); free(lyapunov_avg); fftw_destroy_plan(fftu2.fft_plan); fftw_destroy_plan(fftu3.fft_plan); fftw_destroy_plan(fftu4.fft_plan); fftw_free(fftu2.fft); fftw_free(fftu3.fft); fftw_free(fftu4.fft); ns_free_tmps(u, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9, tmp10, fft1, fftu1, ifft, algorithm); free(tmp11); if(algorithm_lyapunov==ALGORITHM_RK2){ free(tmp1); free(tmp2); } else if (algorithm_lyapunov==ALGORITHM_RK4){ free(tmp1); free(tmp2); free(tmp3); } else { fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm_lyapunov); }; return(0); } // save state of lyapunov computation 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 ){ // save u and step save_state(u, savefile, K1, K2, cmd_string, params_string, savefile_string, utfile_string, utfile, command, algorithm, step, time, nthreads); 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); } return 0; }