From 0f07f025b510c531369183c281fc4695a5ddaa9d Mon Sep 17 00:00:00 2001 From: Ian Jauslin Date: Fri, 31 Jan 2025 10:52:40 -0500 Subject: [PATCH] Overhaul of lyapunov exponents --- src/constants.cpp | 3 + src/lyapunov.c | 726 ++++++++++++++++++++++++++++++++++---------- src/lyapunov.h | 35 ++- src/main.c | 61 +++- src/navier-stokes.c | 3 +- 5 files changed, 651 insertions(+), 177 deletions(-) diff --git a/src/constants.cpp b/src/constants.cpp index 649d916..2eb6422 100644 --- a/src/constants.cpp +++ b/src/constants.cpp @@ -48,3 +48,6 @@ limitations under the License. #define FIX_ENSTROPHY 1 #define FIX_ENERGY 2 + +#define LYAPUNOV_TRIGGER_TIME 1 +#define LYAPUNOV_TRIGGER_SIZE 2 diff --git a/src/lyapunov.c b/src/lyapunov.c index 4d1c3aa..05bcbd5 100644 --- a/src/lyapunov.c +++ b/src/lyapunov.c @@ -32,110 +32,113 @@ int lyapunov( int N2, double final_time, double lyapunov_reset, + unsigned int lyapunov_trigger, double nu, - double D_epsilon, double delta, double L, double adaptive_tolerance, double adaptive_factor, double max_delta, - unsigned int adaptive_norm, + 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* lyapunov; double* lyapunov_avg; - double* Du_prod; - double* Du; + double* flow; _Complex double* u; - _Complex double* prevu; - _Complex double* tmp1; - _Complex double* tmp2; - _Complex double* tmp3; + 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; - double* tmp10; - double time; - fft_vect fft1; - fft_vect fft2; - fft_vect ifft; + _Complex double* tmp10; + double* tmp11; + int i,j; - // period + // 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, &Du_prod, &Du, &u, &prevu, &tmp1, &tmp2, &tmp3, &tmp4, &tmp5, &tmp6, &tmp7, &tmp8, &tmp9, &tmp10, &fft1, &fft2, &ifft, K1, K2, N1, N2, nthreads, algorithm); + 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 Du_prod - int i,j; - for(i=0;i(n+1)*lyapunov_reset){ + double norm=0; + if(lyapunov_trigger==LYAPUNOV_TRIGGER_SIZE){ + // size of flow (for reset) + for(i=0;ilyapunov_reset){ + break; + } + } + if(norm>lyapunov_reset){ + break; + } + } + } + + // QR decomposition + if((lyapunov_trigger==LYAPUNOV_TRIGGER_TIME && time>(n+1)*lyapunov_reset) || (lyapunov_trigger==LYAPUNOV_TRIGGER_SIZE && norm>lyapunov_reset)){ n++; // QR decomposition // do not touch tmp1, it might be used in the next step - LAPACKE_dgerqf(LAPACK_COL_MAJOR, MATSIZE, MATSIZE, Du_prod, MATSIZE, tmp10); - // extract eigenvalues (diagonal elements of Du_prod) + LAPACKE_dgeqrf(LAPACK_COL_MAJOR, MATSIZE, MATSIZE, flow, MATSIZE, tmp11); + // extract diagonal elements of R (represented as diagonal elements of flow for(i=0; i0){ fprintf(stderr," % .15e % .15e\n", lyapunov[0], lyapunov[MATSIZE-1]); } - // set Du_prod to Q: - LAPACKE_dorgrq(LAPACK_COL_MAJOR, MATSIZE, MATSIZE, MATSIZE, Du_prod, MATSIZE, tmp10); + // set flow to Q: + LAPACKE_dorgqr(LAPACK_COL_MAJOR, MATSIZE, MATSIZE, MATSIZE, flow, MATSIZE, tmp11); // reset prevtime prevtime=time; } } - lyapunov_free_tmps(lyapunov, lyapunov_avg, Du_prod, Du, u, prevu, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9, tmp10, fft1, fft2, ifft, algorithm); + 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) { @@ -183,87 +187,64 @@ int compare_double(const void* x, const void* y) { } } -// Jacobian of u_n -> u_{n+1} -int ns_D_step( - double* out, - _Complex double* un, - _Complex double* un_next, +// compute tangent flow +int lyapunov_D_step( + double* flow, + _Complex double* u, 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 fftu1, + fft_vect fftu2, + fft_vect fftu3, + fft_vect fftu4, fft_vect fft1, - fft_vect fft2, fft_vect ifft, - _Complex double* tmp1, - _Complex double* tmp2, - _Complex double* tmp3, + double* tmp1, + double* tmp2, + 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 + bool irreversible ){ int lx,ly; - int i; - double step, next_step; + double alpha; + // 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=-K1;lx<=K1;lx++){ + for(ly=-K2;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++){ - // derivative in the real direction - // finite difference vector - 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)); - *Du_prod=calloc(MATSIZE*MATSIZE,sizeof(double)); - *Du=calloc(MATSIZE*MATSIZE,sizeof(double)); - *prevu=calloc(USIZE,sizeof(_Complex double)); - *tmp1=calloc(USIZE,sizeof(_Complex double)); - *tmp2=calloc(USIZE,sizeof(_Complex double)); - *tmp10=calloc(MATSIZE*MATSIZE,sizeof(double)); + *flow=calloc(MATSIZE*MATSIZE,sizeof(double)); return(0); } // release vectors int lyapunov_free_tmps( - double* lyapunov, - double* lyapunov_avg, - double* Du_prod, - double* Du, - _Complex double* u, - _Complex double* prevu, - _Complex double* tmp1, - _Complex double* tmp2, - _Complex double* tmp3, - _Complex double* tmp4, - _Complex double* tmp5, - _Complex double* tmp6, - _Complex double* tmp7, - _Complex double* tmp8, - _Complex double* tmp9, - double* tmp10, + 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 fft2, fft_vect ifft, - unsigned int algorithm + unsigned int algorithm, + unsigned int algorithm_lyapunov ){ - free(tmp10); - free(tmp2); - free(tmp1); - free(prevu); - free(Du); - free(Du_prod); - free(lyapunov_avg); + free(flow); free(lyapunov); - ns_free_tmps(u, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9, fft1, fft2, ifft, algorithm); + 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); } diff --git a/src/lyapunov.h b/src/lyapunov.h index 9fbbec8..a90c638 100644 --- a/src/lyapunov.h +++ b/src/lyapunov.h @@ -20,18 +20,41 @@ limitations under the License. #include "navier-stokes.h" // 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); +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_norm, _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); // comparison function for qsort int compare_double(const void* x, const void* y); -// Jacobian of step -int ns_D_step( 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); +// 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); + +// 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); + +// RK 2 algorithm +int lyapunov_D_step_rk2( 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, bool irreversible); + +// Right side of equation for tangent flow +int lyapunov_D_rhs( double* out, double* flow, _Complex double* u, int K1, int K2, int N1, int N2, double alpha, 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, bool irreversible); + +// 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); + +// compute derivative of T +int lyapunov_D_T( double* flow, int K1, int K2, int N1, int N2, fft_vect fftu1, fft_vect fftu2, fft_vect fftu3, fft_vect fftu4, fft_vect fft1, fft_vect ifft); + +double lyapunov_D_alpha( double* flow, _Complex double* u, int K1, int K2, int N1, int N2, double alpha, double L, _Complex double* g, _Complex double* T, _Complex double* DT); + +// 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); + +// 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 lyapunov_init_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, int K1, int K2, int N1, int N2, unsigned int nthreads, unsigned int algorithm, unsigned int algorithm_lyapunov); -// init vectors -int lyapunov_init_tmps(double** lyapunov, double** lyapunov_avg, double ** Du_prod, double ** Du, _Complex double ** u, _Complex double ** prevu, _Complex double ** tmp1, _Complex double ** tmp2, _Complex double ** tmp3, _Complex double ** tmp4, _Complex double ** tmp5, _Complex double ** tmp6, _Complex double ** tmp7, _Complex double ** tmp8, _Complex double ** tmp9, double** tmp10, fft_vect* fft1, fft_vect* fft2, fft_vect* ifft, int K1, int K2, int N1, int N2, unsigned int nthreads, unsigned int algorithm); // release vectors -int lyapunov_free_tmps(double* lyapunov, double* lyapunov_avg, double* Du_prod, double* Du, _Complex double* u, _Complex double* prevu, _Complex double* tmp1, _Complex double* tmp2, _Complex double* tmp3, _Complex double* tmp4, _Complex double* tmp5, _Complex double* tmp6, _Complex double* tmp7, _Complex double* tmp8, _Complex double* tmp9, double* tmp10, fft_vect fft1, fft_vect fft2, fft_vect ifft, unsigned int algorithm); +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); #endif diff --git a/src/main.c b/src/main.c index 6ea2216..f7adc08 100644 --- a/src/main.c +++ b/src/main.c @@ -56,10 +56,12 @@ typedef struct nstrophy_parameters { unsigned int driving; unsigned int init; unsigned int algorithm; + unsigned int algorithm_lyapunov; bool keep_en_cst; FILE* initfile; FILE* drivingfile; double lyapunov_reset; + unsigned int lyapunov_trigger; double D_epsilon; bool print_alpha; } nstrophy_parameters; @@ -296,7 +298,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.nu, parameters.D_epsilon, 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.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); } else if(command==0){ fprintf(stderr, "error: no command specified\n"); @@ -576,6 +578,9 @@ int set_default_params( parameters->initfile=NULL; parameters->algorithm=ALGORITHM_RK4; parameters->keep_en_cst=false; + 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->print_alpha=false; @@ -648,6 +653,12 @@ int read_params( parameters->N2=smallest_pow2(3*(parameters->K2)); } + // if lyapunov_reset or lyapunov_maxu are not set explicitly + if(parameters->lyapunov_trigger==0){ + parameters->lyapunov_trigger=LYAPUNOV_TRIGGER_TIME; + parameters->lyapunov_reset=parameters->print_freq; + } + return(0); } @@ -841,13 +852,6 @@ int set_parameter( return(-1); } } - else if (strcmp(lhs,"lyapunov_reset")==0){ - ret=sscanf(rhs,"%lf",&(parameters->lyapunov_reset)); - if(ret!=1){ - fprintf(stderr, "error: parameter 'lyapunov_reset' should be a double\n got '%s'\n",rhs); - return(-1); - } - } else if (strcmp(lhs,"D_epsilon")==0){ ret=sscanf(rhs,"%lf",&(parameters->D_epsilon)); if(ret!=1){ @@ -940,6 +944,47 @@ int set_parameter( } parameters->print_alpha=(tmp==1); } + else if (strcmp(lhs,"lyapunov_reset")==0){ + if(parameters->lyapunov_trigger==LYAPUNOV_TRIGGER_SIZE){ + fprintf(stderr, "error: cannot use 'lyapunov_reset' and 'lyapunov_maxu' in the same run:\n 'lyapunov_reset' is to be used to renormalize the tangent flow at fixed times, 'lyapunov_maxu' is to be used to renormalize the tangent flow when the matrix exceeds a certain size."); + return(-1); + } + parameters->lyapunov_trigger=LYAPUNOV_TRIGGER_TIME; + ret=sscanf(rhs,"%lf",&(parameters->lyapunov_reset)); + if(ret!=1){ + fprintf(stderr, "error: parameter 'lyapunov_reset' should be a double\n got '%s'\n",rhs); + return(-1); + } + } + else if (strcmp(lhs,"lyapunov_maxu")==0){ + if(parameters->lyapunov_trigger==LYAPUNOV_TRIGGER_TIME){ + fprintf(stderr, "error: cannot use 'lyapunov_maxu' and 'lyapunov_reset' in the same run:\n 'lyapunov_reset' is to be used to renormalize the tangent flow at fixed times, 'lyapunov_maxu' is to be used to renormalize the tangent flow when the matrix exceeds a certain size."); + return(-1); + } + parameters->lyapunov_trigger=LYAPUNOV_TRIGGER_SIZE; + ret=sscanf(rhs,"%lf",&(parameters->lyapunov_reset)); + if(ret!=1){ + fprintf(stderr, "error: parameter 'lyapunov_maxu' should be a double\n got '%s'\n",rhs); + return(-1); + } + } + // algorithm for tangent flow (for lyapunov) + else if (strcmp(lhs,"algorithm_lyapunov")==0){ + if (strcmp(rhs,"RK4")==0){ + parameters->algorithm_lyapunov=ALGORITHM_RK4; + } + else if (strcmp(rhs,"RK2")==0){ + parameters->algorithm_lyapunov=ALGORITHM_RK2; + } + else if (strcmp(rhs,"RKF45")==0 || strcmp(rhs,"RKDP45")==0 || strcmp(rhs,"RKBS32")==0){ + fprintf(stderr, "error: cannot use an adaptove step algorithm for the tangent flow (Lyapunov exponents); must be one of 'RK4' or 'RK2', got:'%s'\n",rhs); + return(-1); + } + else{ + fprintf(stderr, "error: unrecognized algorithm '%s'\n",rhs); + return(-1); + } + } else{ fprintf(stderr, "error: unrecognized parameter '%s'\n",lhs); return(-1); diff --git a/src/navier-stokes.c b/src/navier-stokes.c index 873ec59..2f5ee82 100644 --- a/src/navier-stokes.c +++ b/src/navier-stokes.c @@ -82,7 +82,7 @@ int uk( // add 0.1 to ensure proper rounding uint64_t n=(uint64_t)((starting_time-fmod(starting_time, print_freq))/print_freq+0.1); - // store step (useful for adaptive step methods + // store step (useful for adaptive step methods) double step=delta; double next_step=step; @@ -524,6 +524,7 @@ int ns_step( return (0); } +// TODO: do not need to use klookup in any of the rk routines // RK 4 algorithm int ns_step_rk4( _Complex double* u,