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

Split real and imag part in jacobian

Browse Source
This commit is contained in:
Ian Jauslin 2024-03-11 23:54:11 +01:00
parent 9059a24c23
commit 198867751d
3 changed files with 89 additions and 70 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

89
docs/nstrophy_doc.tex
View File

@ -382,62 +382,81 @@ The enstrophy is defined as
To compute the Lyapunov exponents, we must first compute the Jacobian of $\hat u^{(n)}\mapsto\hat u^{(n+1)}$. To compute the Lyapunov exponents, we must first compute the Jacobian of $\hat u^{(n)}\mapsto\hat u^{(n+1)}$.
This map is always of Runge-Kutta type, that is, This map is always of Runge-Kutta type, that is,
\begin{equation} \begin{equation}
\hat u^{(n+1)}=\hat u^{(n)}+\delta\sum_{i=1}^q w_i\mathfrak F(\hat u^{(n)}) \hat u(t_{n+1})=\mathfrak F_{t_n}(\hat u(t_n))
\end{equation}
(see\-~(\ref{ins_k})), so
\begin{equation}
D\hat u^{(n+1)}=\mathds 1+\delta\sum_{i=1}^q w_iD\mathfrak F(\hat u^{(n)})
. .
\end{equation} \end{equation}
We then compute Let $D\mathfrak F_{t_{n}}$ be the Jacobian of this map, in which we split the real and imaginary parts: if
\begin{equation} \begin{equation}
(D\mathfrak F(\hat u))_{k,\ell} \hat u_k(t_n)=:\rho_k+i\iota_k
= ,\quad
-\frac{4\pi^2}{L^2}\nu k^2\delta_{k,\ell} \mathfrak F_{t_n}(\hat u(t_n))_k=:\phi_k+i\psi_k
+\frac{4\pi^2}{L^2|k|}\partial_{\hat u_\ell}T(\hat u,k)
\end{equation} \end{equation}
and, by\-~(\ref{T}), then
\begin{equation} \begin{equation}
\partial_{\hat u_\ell}T(\hat u,k) (D\mathfrak F_{t_n})_{k,p}:=\left(\begin{array}{cc}
= \partial_{\rho_p}\phi_k&\partial_{\iota_p}\phi_k\\
\sum_{\displaystyle\mathop{\scriptstyle q\in\mathbb Z^2}_{\ell+q=k}} \partial_{\rho_p}\psi_k&\partial_{\iota_p}\psi_k
\left( \end{array}\right)
\frac{(q\cdot \ell^\perp)|q|}{|\ell|} \end{equation}
+ We compute this Jacobian numerically using a finite difference, by computing
\frac{(\ell\cdot q^\perp)|\ell|}{|q|} \begin{equation}
\right)\hat u_q (D\mathfrak F_{t_n})_{k,p}:=\frac1\epsilon
= \left(\begin{array}{cc}
(k\cdot \ell^\perp)\left( \phi_k(\hat u+\epsilon\delta_p)-\phi_k(\hat u)&\phi_k(\hat u+i\epsilon\delta_p)-\phi_k(\hat u)\\
\frac{|k-\ell|}{|\ell|} \psi_k(\hat u+\epsilon\delta_p)-\psi_k(\hat u)&\psi_k(\hat u+i\epsilon\delta_p)-\psi_k(\hat u)
- \end{array}\right)
\frac{|\ell|}{|k-\ell|}
\right)\hat u_{k-\ell}
. .
\end{equation} \end{equation}
The parameter $\epsilon$ can be set using the parameter {\tt D\_epsilon}.
%, so
%\begin{equation}
% D\hat u^{(n+1)}=\mathds 1+\delta\sum_{i=1}^q w_iD\mathfrak F(\hat u^{(n)})
% .
%\end{equation}
%We then compute
%\begin{equation}
% (D\mathfrak F(\hat u))_{k,\ell}
% =
% -\frac{4\pi^2}{L^2}\nu k^2\delta_{k,\ell}
% +\frac{4\pi^2}{L^2|k|}\partial_{\hat u_\ell}T(\hat u,k)
%\end{equation}
%and, by\-~(\ref{T}),
%\begin{equation}
% \partial_{\hat u_\ell}T(\hat u,k)
% =
% \sum_{\displaystyle\mathop{\scriptstyle q\in\mathbb Z^2}_{\ell+q=k}}
% \left(
% \frac{(q\cdot \ell^\perp)|q|}{|\ell|}
% +
% \frac{(\ell\cdot q^\perp)|\ell|}{|q|}
% \right)\hat u_q
% =
% (k\cdot \ell^\perp)\left(
% \frac{|k-\ell|}{|\ell|}
% -
% \frac{|\ell|}{|k-\ell|}
% \right)\hat u_{k-\ell}
% .
%\end{equation}
\bigskip \bigskip
\indent \indent
The Lyapunov exponents are then defined as The Lyapunov exponents are then defined as
\begin{equation} \begin{equation}
\frac1n\log\mathrm{spec}(\pi_i) \frac1{t_n}\log\mathrm{spec}(\pi_{t_n})
,\quad ,\quad
\pi_i:=\prod_{i=1}^nD\hat u^{(i)} \pi_{t_n}:=\prod_{i=1}^nD\hat u(t_i)
. .
\end{equation} \end{equation}
However, the product of $D\hat u$ may become quite large or quite small (if the exponents are not all 1). However, the product of $D\hat u$ may become quite large or quite small (if the exponents are not all 1).
To avoid this, we periodically rescale the product. To avoid this, we periodically rescale the product.
We set $\mathfrak L_r\in\mathbb N_*$ (set by adjusting the {\tt lyanpunov\_reset} parameter), and, when $n$ is a multiple of $\mathfrak L_r$, we rescale the eigenvalues of $\pi_i$ to 1. We set $\mathfrak L_r>0$ (set by adjusting the {\tt lyanpunov\_reset} parameter), and, when $t_n$ crosses a multiple of $\mathfrak L_r$, we rescale the eigenvalues of $\pi_i$ to 1.
To do so, we perform a $QR$ decomposition: To do so, we perform a $QR$ decomposition:
\begin{equation} \begin{equation}
\pi_{\alpha\mathfrak L_r} \pi_{\alpha\mathfrak L_r}
=Q^{(\alpha)}R^{(\alpha)} =R^{(\alpha)}Q^{(\alpha)}
\end{equation}
where $Q^{(\alpha)}$ is diagonal and $R^{(\alpha)}$ is an orthogonal matrix, and we divide by $Q^{(\alpha)}$ (thus only keeping $R^{(\alpha)}$).
Thus, we replace
\begin{equation}
\pi_{\alpha\mathfrak L_r+i}\mapsto R^{(\alpha)}\prod_{j=1}^iD\hat u^{(j)}
.
\end{equation} \end{equation}
where $Q^{(\alpha)}$ is orthogonal and $R^{(\alpha)}$ is a diagonal matrix, and we divide by $R^{(\alpha)}$ (thus only keeping $Q^{(\alpha)}$).
The Lyapunov exponents at time $\alpha\mathfrak L_r$ are then The Lyapunov exponents at time $\alpha\mathfrak L_r$ are then
\begin{equation} \begin{equation}
\frac1{\alpha\mathfrak L_r}\sum_{\beta=1}^\alpha\log\mathrm{spec}(Q^{(\beta)}) \frac1{\alpha\mathfrak L_r}\sum_{\beta=1}^\alpha\log\mathrm{spec}(Q^{(\beta)})

64
src/lyapunov.c
View File

@ -21,7 +21,8 @@ limitations under the License.
#include <math.h> #include <math.h>
#include <stdlib.h> #include <stdlib.h>
#define MATSIZE (K1*(2*K2+1)+K2) #define MATSIZE (2*(K1*(2*K2+1)+K2))
#define USIZE (K1*(2*K2+1)+K2)
// compute Lyapunov exponents // compute Lyapunov exponents
int lyapunov( int lyapunov(
@ -50,8 +51,8 @@ int lyapunov(
){ ){
double* lyapunov; double* lyapunov;
double* lyapunov_avg; double* lyapunov_avg;
_Complex double* Du_prod; double* Du_prod;
_Complex double* Du; double* Du;
_Complex double* u; _Complex double* u;
_Complex double* prevu; _Complex double* prevu;
_Complex double* tmp1; _Complex double* tmp1;
@ -63,7 +64,7 @@ int lyapunov(
_Complex double* tmp7; _Complex double* tmp7;
_Complex double* tmp8; _Complex double* tmp8;
_Complex double* tmp9; _Complex double* tmp9;
_Complex double* tmp10; double* tmp10;
double time; double time;
fft_vect fft1; fft_vect fft1;
fft_vect fft2; fft_vect fft2;
@ -82,9 +83,6 @@ int lyapunov(
double step=delta; double step=delta;
double next_step=step; double next_step=step;
const _Complex double one=1;
const _Complex double zero=0;
int i; int i;
// init average // init average
for (i=0; i<MATSIZE; i++){ for (i=0; i<MATSIZE; i++){
@ -109,9 +107,9 @@ int lyapunov(
// multiply Jacobian // multiply Jacobian
// switch to column major order, so transpose Du // 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); cblas_dgemm(CblasColMajor, CblasNoTrans, CblasTrans, MATSIZE, MATSIZE, MATSIZE, 1., Du_prod, MATSIZE, Du, MATSIZE, 0., tmp10, MATSIZE);
// copy back to Du_prod // copy back to Du_prod
_Complex double* move=tmp10; double* move=tmp10;
tmp10=Du_prod; tmp10=Du_prod;
Du_prod=move; Du_prod=move;
@ -124,10 +122,10 @@ int lyapunov(
// QR decomposition // QR decomposition
// do not touch tmp1, it might be used in the next step // do not touch tmp1, it might be used in the next step
LAPACKE_zgerqf(LAPACK_COL_MAJOR, MATSIZE, MATSIZE, Du_prod, MATSIZE, tmp2); LAPACKE_dgerqf(LAPACK_COL_MAJOR, MATSIZE, MATSIZE, Du_prod, MATSIZE, tmp10);
// extract eigenvalues (diagonal elements of Du_prod) // extract eigenvalues (diagonal elements of Du_prod)
for(i=0; i<MATSIZE; i++){ for(i=0; i<MATSIZE; i++){
lyapunov[i]=log(cabs(Du_prod[i*MATSIZE+i]))/(time-prevtime); lyapunov[i]=log(fabs(Du_prod[i*MATSIZE+i]))/(time-prevtime);
} }
// sort lyapunov exponents // sort lyapunov exponents
@ -153,7 +151,7 @@ int lyapunov(
} }
// set Du_prod to Q: // set Du_prod to Q:
LAPACKE_zungrq(LAPACK_COL_MAJOR, MATSIZE, MATSIZE, MATSIZE, Du_prod, MATSIZE, tmp2); LAPACKE_dorgrq(LAPACK_COL_MAJOR, MATSIZE, MATSIZE, MATSIZE, Du_prod, MATSIZE, tmp10);
// reset prevtime // reset prevtime
prevtime=time; prevtime=time;
@ -177,7 +175,7 @@ int compare_double(const void* x, const void* y) {
// Jacobian of u_n -> u_{n+1} // Jacobian of u_n -> u_{n+1}
int ns_D_step( int ns_D_step(
_Complex double* out, double* out,
_Complex double* un, _Complex double* un,
_Complex double* un_next, _Complex double* un_next,
int K1, int K1,
@ -218,7 +216,7 @@ int ns_D_step(
for(ly=(lx>0 ? -K2 : 1);ly<=K2;ly++){ for(ly=(lx>0 ? -K2 : 1);ly<=K2;ly++){
// derivative in the real direction // derivative in the real direction
// finite difference vector // finite difference vector
for (i=0;i<MATSIZE;i++){ for (i=0;i<USIZE;i++){
if(i==klookup_sym(lx,ly,K2)){ if(i==klookup_sym(lx,ly,K2)){
tmp1[i]=un[i]+epsilon; tmp1[i]=un[i]+epsilon;
}else{ }else{
@ -230,15 +228,16 @@ int ns_D_step(
step=delta; step=delta;
next_step=delta; next_step=delta;
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); 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 // derivative
for (i=0;i<MATSIZE;i++){ for (i=0;i<USIZE;i++){
// use row major order // use row major order
out[klookup_sym(lx,ly,K2)*MATSIZE+i]=(tmp1[i]-un_next[i])/epsilon; out[2*klookup_sym(lx,ly,K2)*USIZE+2*i]=(__real__ (tmp1[i]-un_next[i]))/epsilon;
out[(2*klookup_sym(lx,ly,K2)+1)*USIZE+2*i]=(__imag__ (tmp1[i]-un_next[i]))/epsilon;
} }
// derivative in the imaginary direction // derivative in the imaginary direction
// finite difference vector // finite difference vector
for (i=0;i<MATSIZE;i++){ for (i=0;i<USIZE;i++){
if(i==klookup_sym(lx,ly,K2)){ if(i==klookup_sym(lx,ly,K2)){
tmp1[i]=un[i]+epsilon*I; tmp1[i]=un[i]+epsilon*I;
}else{ }else{
@ -251,9 +250,10 @@ int ns_D_step(
next_step=delta; next_step=delta;
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); 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 // compute derivative
for (i=0;i<MATSIZE;i++){ for (i=0;i<USIZE;i++){
// use row major order // use row major order
out[klookup_sym(lx,ly,K2)*MATSIZE+i]=-(tmp1[i]-un_next[i])/epsilon*I; out[2*klookup_sym(lx,ly,K2)*USIZE+2*i+1]=(__real__ (tmp1[i]-un_next[i]))/epsilon;
out[(2*klookup_sym(lx,ly,K2)+1)*USIZE+2*i+1]=(__imag__ (tmp1[i]-un_next[i]))/epsilon;
} }
} }
} }
@ -265,8 +265,8 @@ int ns_D_step(
int lyapunov_init_tmps( int lyapunov_init_tmps(
double ** lyapunov, double ** lyapunov,
double ** lyapunov_avg, double ** lyapunov_avg,
_Complex double ** Du_prod, double ** Du_prod,
_Complex double ** Du, double ** Du,
_Complex double ** u, _Complex double ** u,
_Complex double ** prevu, _Complex double ** prevu,
_Complex double ** tmp1, _Complex double ** tmp1,
@ -278,7 +278,7 @@ int lyapunov_init_tmps(
_Complex double ** tmp7, _Complex double ** tmp7,
_Complex double ** tmp8, _Complex double ** tmp8,
_Complex double ** tmp9, _Complex double ** tmp9,
_Complex double ** tmp10, double ** tmp10,
fft_vect* fft1, fft_vect* fft1,
fft_vect* fft2, fft_vect* fft2,
fft_vect* ifft, fft_vect* ifft,
@ -293,12 +293,12 @@ int lyapunov_init_tmps(
*lyapunov=calloc(sizeof(double),MATSIZE); *lyapunov=calloc(sizeof(double),MATSIZE);
*lyapunov_avg=calloc(sizeof(double),MATSIZE); *lyapunov_avg=calloc(sizeof(double),MATSIZE);
*Du_prod=calloc(sizeof(_Complex double),MATSIZE*MATSIZE); *Du_prod=calloc(sizeof(double),MATSIZE*MATSIZE);
*Du=calloc(sizeof(_Complex double),MATSIZE*MATSIZE); *Du=calloc(sizeof(double),MATSIZE*MATSIZE);
*prevu=calloc(sizeof(_Complex double),MATSIZE); *prevu=calloc(sizeof(_Complex double),USIZE);
*tmp1=calloc(sizeof(_Complex double),MATSIZE); *tmp1=calloc(sizeof(_Complex double),USIZE);
*tmp2=calloc(sizeof(_Complex double),MATSIZE); *tmp2=calloc(sizeof(_Complex double),USIZE);
*tmp10=calloc(sizeof(_Complex double),MATSIZE*MATSIZE); *tmp10=calloc(sizeof(double),MATSIZE*MATSIZE);
return(0); return(0);
} }
@ -307,8 +307,8 @@ int lyapunov_init_tmps(
int lyapunov_free_tmps( int lyapunov_free_tmps(
double* lyapunov, double* lyapunov,
double* lyapunov_avg, double* lyapunov_avg,
_Complex double* Du_prod, double* Du_prod,
_Complex double* Du, double* Du,
_Complex double* u, _Complex double* u,
_Complex double* prevu, _Complex double* prevu,
_Complex double* tmp1, _Complex double* tmp1,
@ -320,7 +320,7 @@ int lyapunov_free_tmps(
_Complex double* tmp7, _Complex double* tmp7,
_Complex double* tmp8, _Complex double* tmp8,
_Complex double* tmp9, _Complex double* tmp9,
_Complex double* tmp10, double* tmp10,
fft_vect fft1, fft_vect fft1,
fft_vect fft2, fft_vect fft2,
fft_vect ifft, fft_vect ifft,

6
src/lyapunov.h
View File

@ -26,12 +26,12 @@ int lyapunov( int K1, int K2, int N1, int N2, double final_time, double lyapunov
int compare_double(const void* x, const void* y); int compare_double(const void* x, const void* y);
// Jacobian of step // 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); 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);
// init vectors // 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); 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 // 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); 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);
#endif #endif