Computation of Lyapunov exponents
This commit is contained in:
parent
89601d19d1
commit
b31fab8eae
208
src/lyapunov.c
208
src/lyapunov.c
@ -16,10 +16,134 @@ limitations under the License.
|
||||
|
||||
#include "constants.cpp"
|
||||
#include "lyapunov.h"
|
||||
#include <cblas.h>
|
||||
#include <openblas64/cblas.h>
|
||||
#include <openblas64/lapacke.h>
|
||||
#include <math.h>
|
||||
#include <stdlib.h>
|
||||
|
||||
#define MATSIZE (K1*(2*(K2+1)+K2))
|
||||
|
||||
// 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
|
||||
){
|
||||
double* lyapunov;
|
||||
double* lyapunov_avg;
|
||||
_Complex double* Du_prod;
|
||||
_Complex 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;
|
||||
_Complex double* tmp10;
|
||||
double time;
|
||||
fft_vect fft1;
|
||||
fft_vect fft2;
|
||||
fft_vect ifft;
|
||||
|
||||
// period
|
||||
// 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);
|
||||
// copy initial condition
|
||||
copy_u(u, u0, K1, K2);
|
||||
|
||||
// store step (useful for adaptive step methods
|
||||
double prev_step=delta;
|
||||
double step=delta;
|
||||
double next_step=step;
|
||||
|
||||
const _Complex double one=1;
|
||||
const _Complex double zero=0;
|
||||
|
||||
int i;
|
||||
// init average
|
||||
for (i=0; i<MATSIZE; i++){
|
||||
lyapunov_avg=0;
|
||||
}
|
||||
|
||||
// save times at which Lyapunov exponents are computed
|
||||
double prevtime=starting_time;
|
||||
|
||||
// iterate
|
||||
time=starting_time;
|
||||
while(final_time<0 || time<final_time){
|
||||
// copy before step
|
||||
copy_u(prevu, u, K1, K2);
|
||||
prev_step=step;
|
||||
|
||||
ns_step(algorithm, u, K1, K2, N1, N2, nu, &step, &next_step, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, L, g, time, starting_time, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en);
|
||||
|
||||
// compute Jacobian
|
||||
// do not touch tmp1, it might be used in the next step
|
||||
ns_D_step(Du, prevu, u, K1, K2, N1, N2, nu, D_epsilon, prev_step, algorithm, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, L, g, time, fft1, fft2, ifft, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9, irreversible, keep_en_cst, target_en);
|
||||
|
||||
// multiply Jacobian
|
||||
// switch to column major order, so transpose Du
|
||||
cblas_zgemm(CblasColMajor, CblasNoTrans, CblasTrans, MATSIZE, MATSIZE, MATSIZE, &one, Du_prod, MATSIZE, Du, MATSIZE, &zero, tmp10, MATSIZE);
|
||||
// copy back to Du_prod
|
||||
_Complex double* move=tmp10;
|
||||
tmp10=Du_prod;
|
||||
Du_prod=move;
|
||||
|
||||
// increment time
|
||||
time+=step;
|
||||
|
||||
// reset Jacobian
|
||||
if(time>(n+1)*lyapunov_reset){
|
||||
n++;
|
||||
|
||||
// QR decomposition
|
||||
// do not touch tmp1, it might be used in the next step
|
||||
LAPACKE_zgerqf(LAPACK_COL_MAJOR, MATSIZE, MATSIZE, Du_prod, MATSIZE, tmp2);
|
||||
// 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
|
||||
lyapunov_avg[i]=lyapunov_avg[i]*prevtime/time+lyapunov[i]*(time-prevtime)/time;
|
||||
}
|
||||
|
||||
// set Du_prod to Q:
|
||||
LAPACKE_zungrq(LAPACK_COL_MAJOR, MATSIZE, MATSIZE, MATSIZE, Du_prod, MATSIZE, tmp2);
|
||||
|
||||
// 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);
|
||||
return(0);
|
||||
}
|
||||
|
||||
// Jacobian of u_n -> u_{n+1}
|
||||
int ns_D_step(
|
||||
_Complex double* out,
|
||||
@ -77,7 +201,8 @@ int ns_D_step(
|
||||
ns_step(algorithm, tmp1, K1, K2, N1, N2, nu, &step, &next_step, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, L, g, time, time, fft1, fft2, ifft, &tmp2, &tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, irreversible, keep_en_cst, target_en);
|
||||
// compute derivative
|
||||
for (i=0;i<MATSIZE;i++){
|
||||
out[i]=(tmp1[i]-un_next[i])/epsilon;
|
||||
// use row major order
|
||||
out[klookup_sym(lx,ly,K2)*MATSIZE+i]=(tmp1[i]-un_next[i])/epsilon;
|
||||
}
|
||||
|
||||
// derivative in the imaginary direction
|
||||
@ -96,7 +221,8 @@ int ns_D_step(
|
||||
ns_step(algorithm, tmp1, K1, K2, N1, N2, nu, &step, &next_step, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, L, g, time, time, fft1, fft2, ifft, &tmp2, &tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, irreversible, keep_en_cst, target_en);
|
||||
// compute derivative
|
||||
for (i=0;i<MATSIZE;i++){
|
||||
out[i]+=-(tmp1[i]-un_next[i])/epsilon*I;
|
||||
// use row major order
|
||||
out[klookup_sym(lx,ly,K2)*MATSIZE+i]=-(tmp1[i]-un_next[i])/epsilon*I;
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -105,3 +231,79 @@ int ns_D_step(
|
||||
}
|
||||
|
||||
|
||||
int lyapunov_init_tmps(
|
||||
double ** lyapunov,
|
||||
double ** lyapunov_avg,
|
||||
_Complex double ** Du_prod,
|
||||
_Complex 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,
|
||||
_Complex 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
|
||||
){
|
||||
ns_init_tmps(u, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9, fft1, fft2, ifft, K1, K2, N1, N2, nthreads, algorithm);
|
||||
|
||||
*lyapunov=calloc(sizeof(double),MATSIZE);
|
||||
*lyapunov_avg=calloc(sizeof(double),MATSIZE);
|
||||
*Du_prod=calloc(sizeof(_Complex double),MATSIZE*MATSIZE);
|
||||
*Du=calloc(sizeof(_Complex double),MATSIZE*MATSIZE);
|
||||
*prevu=calloc(sizeof(_Complex double),MATSIZE);
|
||||
*tmp1=calloc(sizeof(_Complex double),MATSIZE);
|
||||
*tmp2=calloc(sizeof(_Complex double),MATSIZE);
|
||||
*tmp10=calloc(sizeof(_Complex double),MATSIZE*MATSIZE);
|
||||
|
||||
return(0);
|
||||
}
|
||||
|
||||
// release vectors
|
||||
int lyapunov_free_tmps(
|
||||
double* lyapunov,
|
||||
double* lyapunov_avg,
|
||||
_Complex double* Du_prod,
|
||||
_Complex 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,
|
||||
_Complex double* tmp10,
|
||||
fft_vect fft1,
|
||||
fft_vect fft2,
|
||||
fft_vect ifft,
|
||||
unsigned int algorithm
|
||||
){
|
||||
free(tmp10);
|
||||
free(tmp2);
|
||||
free(tmp1);
|
||||
free(prevu);
|
||||
free(Du);
|
||||
free(Du_prod);
|
||||
free(lyapunov_avg);
|
||||
free(lyapunov);
|
||||
ns_free_tmps(u, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9, fft1, fft2, ifft, algorithm);
|
||||
|
||||
return(0);
|
||||
}
|
||||
|
||||
@ -19,7 +19,17 @@ 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);
|
||||
|
||||
// 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);
|
||||
|
||||
// init vectors
|
||||
int lyapunov_init_tmps(double** lyapunov, double** lyapunov_avg, _Complex double ** Du_prod, _Complex 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, _Complex 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, _Complex double* Du_prod, _Complex 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, _Complex double* tmp10, fft_vect fft1, fft_vect fft2, fft_vect ifft, unsigned int algorithm);
|
||||
|
||||
#endif
|
||||
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user