Remove init_enstrophy and init_energy from savefile

Add command to print the enstrophy of the initial condition
Fix remove_entry
2025-12-14 10:29:47 -05:00 · 2025-10-14 12:20:43 -04:00 · 2025-06-05 14:55:34 +02:00 · 2025-04-14 18:25:59 -04:00 · 2025-03-04 10:57:05 -05:00 · 2025-02-10 15:34:03 -05:00
21 changed files with 1425 additions and 549 deletions
--- a/2
+++ b/2
@@ -1,4 +1,4 @@
-# Copyright 2017-2024 Ian Jauslin
+# 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.
--- a/README.md
+++ b/README.md
@@ -39,15 +39,33 @@ The available commands are

 * `enstrophy`: to compute the enstrophy and various other observables. This
  command prints
-  ```step_index time average(alpha) average(energy) average(enstrophy) alpha energy enstrophy```
+  ```
+  step_index time average(alpha) average(energy) average(enstrophy) alpha energy enstrophy Re(u(1,1)) Re(u(1,2))
+  ```
  where the averages are running averages over `print_freq`. In addition, if
  the algorithm has an adaptive step, an extra column is printed with `delta`.
  In addition, if alpha has a negative value (even in between `print_freq`
-  intervals), a line is printed with the information.
+  intervals), a line is printed with the information. The two components (1,1)
+  and (1,2) of u are included to more easily identify periodic trajectories, or
+  the presence of multiple attractors.
+
+* `lyapunov`: to compute the Lyapunov exponents. This command prints
+  ```
+  time instantaneous_lyapunov lyapunov
+  ```
+  where `instantaneous_lyapunov` is computed from the tangent flow only between
+  the given time and the previous one.

 * `uk`: to compute the Fourier transform of the solution.

-* `quiet`: does not print anything, useful for debugging.
+* `quiet`: does not print anything (useful for debugging).
+
+* `enstrophy_print_init`: to compute the enstrophy and various other
+  observables for the initial condition (useful for debugging). The command
+  prints
+  ```
+  alpha energy enstrophy
+  ```


 # Parameters
@@ -138,20 +156,35 @@ should be a `;` sperated list of `key=value` pairs. The possible keys are
 * `print_alpha` (0 or 1, default 0): if this is set to 1, then whenever alpha
  is negative, its value is printed as a comment.

+* `lyapunov_reset` (double, default: `print_freq`): if this is set, then, when
+  computing the Lyapunov exponents, the tangent flow will renormalize itself at
+  times proportional to `lyapunov_reset`. This option is incompatible with
+  `lyapunov_maxu`.
+
+* `lyapunov_maxu` (double, default: unset): if this is set, then, when
+  computing the Lyapunov exponents, the tangent flow will renormalize itself
+  whenever the norm of the tangent flow exceeds  `lyapunov_maxu`.  This option
+  is incompatible with `lyapunov_reset`.
+
+* `algorithm_lyapunov`: the algorithm used to integrate the tangent flow. Can
+  be `RK4` for Runge-Kutta 4 (default) or `RK2` for Runge-Kutta 2. Adaptive
+  step algorithms cannot be used for the tangent flow.
+

 # Interrupting and resuming the computation

-The computation can be interrupted by sending Nstrophy the `SIGINT` signal
-(e.g. by pressing `Ctrl-C`.) When Nstrophy receives the `SIGINT` signal, it
-finishes its current step and writes the value of uk, either to `savefile` if
-such a file was specified on the command line (using the `-s` flag), or to
-`stderr`. In addition, when a `savefile` is specified it writes the command
-that needs to be used to resume the computation (which essentially just sets
-the appropriate `starting_time` and `init:file:<savefile>` parameters. The data
-written to the `savefile` is binary.
+The `enstrophy` and `lyapunov` computations can be interrupted by sending
+Nstrophy the `SIGINT` signal (e.g. by pressing `Ctrl-C`.) When Nstrophy
+receives the `SIGINT` signal, it finishes its current step and writes the value
+of uk, either to `savefile` if such a file was specified on the command line
+(using the `-s` flag), or to `stderr`. In addition, when a `savefile` is
+specified it writes the command that needs to be used to resume the computation
+(which essentially just sets the appropriate `starting_time` and
+`init:file:<savefile>` parameters). The data written to the `savefile` is
+binary.


 # License
 Nstrophy is released under the Apache 2.0 license.

-Copyright 2017-2024 Ian Jauslin
+Copyright 2017-2025 Ian Jauslin
--- a/docs/nstrophy_doc.tex
+++ b/docs/nstrophy_doc.tex
@@ -1,4 +1,4 @@
-% Copyright 2017-2024 Ian Jauslin
+% 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.
@@ -103,9 +103,97 @@ We truncate the Fourier modes and assume that $\hat u_k=0$ if $|k_1|>K_1$ or $|k
  \mathcal K:=\{(k_1,k_2),\ |k_1|\leqslant K_1,\ |k_2|\leqslant K_2\}
  .
 \end{equation}
+\bigskip

-\subsubsection{Runge-Kutta methods}.
-To solve the equation numerically, we will use Runge-Kutta methods, which compute an approximate value $\hat u_k^{(n)}$ for $\hat u_k(t_n)$.
+{\bf Remark}:
+Since $U$ is real, $\hat U_{-k}=\hat U_k^*$, and so
+\begin{equation}
+  \hat u_{-k}=\hat u_k^*
+  .
+  \label{realu}
+\end{equation}
+Similarly,
+\begin{equation}
+  \hat g_{-k}=\hat g_k^*
+  .
+  \label{realg}
+\end{equation}
+Thus,
+\begin{equation}
+  T(\hat u,-k)
+  =
+  T(\hat u,k)^*
+  .
+  \label{realT}
+\end{equation}
+In order to keep the computation as quick as possible, we only compute and store the values for $k_1\geqslant 0$.
+(In fact, if we do not enforce the reality conditions, the computation has been found to be unstable.)
+
+\subsection{Reversible equation}
+\indent The reversible equation is similar to\-~(\ref{ins}) but instead of fixing the viscosity, we fix the enstrophy\-~\cite{Ga22}.
+It is defined directly in Fourier space:
+\begin{equation}
+  \partial_t\hat U_k=
+  -\frac{4\pi^2}{L^2}\alpha(\hat U) k^2\hat U_k+\hat G_k
+  -i\frac{2\pi}L\sum_{\displaystyle\mathop{\scriptstyle p,q\in\mathbb Z^2}_{p+q=k}}
+  (q\cdot\hat U_p)\hat U_q
+  ,\quad
+  k\cdot\hat U_k=0
+\end{equation}
+where $\alpha$ is chosen such that the enstrophy is constant.
+In terms of $\hat u$\-~(\ref{udef}), (\ref{gdef}), (\ref{T}):
+\begin{equation}
+  \partial_t\hat u_k=
+  -\frac{4\pi^2}{L^2}\alpha(\hat u) k^2\hat u_k
+  +\hat g_k
+  +\frac{4\pi^2}{L^2|k|}T(\hat u,k)
+  .
+  \label{rns_k}
+\end{equation}
+To compute $\alpha$, we use the constancy of the enstrophy:
+\begin{equation}
+  \sum_{k\in\mathbb Z^2}k^2\hat U_k\cdot\partial_t\hat U_k
+  =0
+\end{equation}
+which, in terms of $\hat u$ is
+\begin{equation}
+  \sum_{k\in\mathbb Z^2}k^2\hat u_k^*\partial_t\hat u_k
+  =0
+\end{equation}
+that is
+\begin{equation}
+  \frac{4\pi^2}{L^2}\alpha(\hat u)\sum_{k\in\mathbb Z^2}k^4|\hat u_k|^2
+  =
+  \sum_{k\in\mathbb Z^2}k^2\hat u_k^*\hat g_k
+  +\frac{4\pi^2}{L^2}\sum_{k\in\mathbb Z^2}|k|\hat u_k^*T(\hat u,k)
+\end{equation}
+and so
+\begin{equation}
+  \alpha(\hat u)
+  =\frac{\frac{L^2}{4\pi^2}\sum_k k^2\hat u_k^*\hat g_k+\sum_k|k|\hat u_k^*T(\hat u,k)}{\sum_kk^4|\hat u_k|^2}
+  .
+  \label{alpha}
+\end{equation}
+Note that, by\-~(\ref{realu})-(\ref{realT}),
+\begin{equation}
+  \alpha(\hat u)\in\mathbb R
+  .
+\end{equation}
+
+\subsection{Runge-Kutta methods}.
+To solve these equations numerically, we will use Runge-Kutta methods, which compute an approximate value $\hat u_k^{(n)}$ for $\hat u_k(t_n)$.
+These algorithms approximate the solution to an equation of the form $\dot u=f(t;u)$ with
+\begin{equation}
+  \hat u_k^{(n+1)}=\hat u_k^{(n)}
+  +\delta_n\sum_{i=1}^sb_ik_i(t_n;\hat u^{(n)})
+  ,\quad
+  k_i(t_n;\hat u^{(n)}):=f\left(t_n+c_i\delta_n;\ \hat u+\delta_n\sum_{j=1}^{i-1}a_{i,j}k_j(t_n,\hat u^{(n)})\right)
+  .
+\end{equation}
+The $c_i$ and $a_{i,j}$ are chosen in one of various ways, depending on the desired accuracy.
+\bigskip
+
+\indent
 {\tt nstrophy} supports the 4th order Runge-Kutta ({\tt RK4}) and 2nd order Runge-Kutta ({\tt RK2}) algorithms.
 In addition, several variable step methods are implemented:
 \begin{itemize}
@@ -190,30 +278,8 @@ It can be made by specifying the parameter {\tt adaptive\_cost}.
  These cost functions are selected by choosing {\tt adaptive\_cost=k3} and {\tt adaptive\_cost=k32} respectively.
 \end{itemize}

-\subsubsection{Reality}.
-Since $U$ is real, $\hat U_{-k}=\hat U_k^*$, and so
-\begin{equation}
-  \hat u_{-k}=\hat u_k^*
-  .
-  \label{realu}
-\end{equation}
-Similarly,
-\begin{equation}
-  \hat g_{-k}=\hat g_k^*
-  .
-  \label{realg}
-\end{equation}
-Thus,
-\begin{equation}
-  T(\hat u,-k)
-  =
-  T(\hat u,k)^*
-  .
-  \label{realT}
-\end{equation}
-In order to keep the computation as quick as possible, we only compute and store the values for $k_1\geqslant 0$.

-\subsubsection{FFT}. We compute T using a fast Fourier transform, defined as
+\subsection{Computation of $T$: FFT}. We compute T using a fast Fourier transform, defined as
 \begin{equation}
  \mathcal F(f)(n):=\sum_{m\in\mathcal N}e^{-\frac{2i\pi}{N_1}m_1n_1-\frac{2i\pi}{N_2}m_2n_2}f(m_1,m_2)
 \end{equation}
@@ -235,7 +301,7 @@ in which $\mathcal F^*$ is defined like $\mathcal F$ but with the opposite phase
  \sum_{p,q\in\mathcal K}
  \frac1{N_1N_2}
  \sum_{n\in\mathcal N}e^{-\frac{2i\pi}{N_1}n_1(p_1+q_1-k_1)-\frac{2i\pi}{N_2}n_2(p_2+q_2-k_2)}
-  (q\cdot p^\perp)\frac{|q|}{|p|}\hat u_q\hat u_p
+  (q\cdot p^\perp)\frac{|q|}{|p|}\hat u_p\hat u_q
 \end{equation}
 provided
 \begin{equation}
@@ -257,6 +323,11 @@ Therefore,
  \right)(k)
 \end{equation}

+\subsection{Observables}
+\indent
+We define the following observables.
+\bigskip
+
 \subsubsection{Energy}.
 We define the energy as
 \begin{equation}
@@ -367,49 +438,125 @@ The enstrophy is defined as
  .
 \end{equation}

-\subsubsection{Lyapunov exponents}
+\subsection{Lyapunov exponents}
 \indent
-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,
+The Lyapunov are defined from the {\it tangent flow} of the dynamics.
+Consider an equation of the form
 \begin{equation}
-  \hat u(t_{n+1})=\mathfrak F_{t_n}(\hat u(t_n))
+  \dot u=f(t;u)
  .
 \end{equation}
-Let $D\mathfrak F_{t_{n}}$ be the Jacobian of this map, in which we split the real and imaginary parts: if
+Now, the flow may not be complex-differentiable, so the tangent flow should be computed on the real and imaginary parts.
+Let
 \begin{equation}
-  \hat u_k(t_n)=:\rho_k+i\iota_k
+  u=\zeta+i\xi
  ,\quad
-  \mathfrak F_{t_n}(\hat u(t_n))_k=:\phi_k+i\psi_k
-\end{equation}
-then
-\begin{equation}
-  (D\mathfrak F_{t_n})_{k,p}:=\left(\begin{array}{cc}
-    \partial_{\rho_p}\phi_k&\partial_{\iota_p}\phi_k\\
-    \partial_{\rho_p}\psi_k&\partial_{\iota_p}\psi_k
-  \end{array}\right)
-\end{equation}
-We compute this Jacobian numerically using a finite difference, by computing
-\begin{equation}
-  (D\mathfrak F_{t_n})_{k,p}:=\frac1\epsilon
-  \left(\begin{array}{cc}
-    \phi_k(\hat u+\epsilon\delta_p)-\phi_k(\hat u)&\phi_k(\hat u+i\epsilon\delta_p)-\phi_k(\hat u)\\
-    \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)
+  f(t;u)=\theta(t;\zeta,\xi)+i\psi(t;\zeta,\xi)
+  .
+\end{equation}
+The tangent flow is given by
+\begin{equation}
+  \dot\delta=\left(\begin{array}{cc}
+    D_\zeta\theta&D_\xi\theta\\
+    D_\zeta\psi&D_\xi\psi
+  \end{array}\right)\delta
+\end{equation}
+where $D_\zeta\theta$ is the Jacobian of $\theta$ with respect to $\zeta$ and so forth...
+The flow of this equation is denoted by
+\begin{equation}
+  \varphi_{t_0,t_1}(\delta_0)
+\end{equation}
+and defined by
+\begin{equation}
+  \frac d{dt}\varphi_{t_0,t}(\delta_0)=\left(\begin{array}{cc}
+    D_\zeta\theta(t;\zeta,\xi)&D_\xi\theta(t;\zeta,\xi)\\
+    D_\zeta\psi(t;\zeta,\xi)&D_\xi\psi(t;\zeta,\xi)
+  \end{array}\right)\varphi_{t_0,t}(\delta_0)
+  ,\quad
+  \varphi_{t_0,t_0}(\delta_0)=\delta_0
+  .
+\end{equation}
+The flow $\varphi_{t_0,t}(\delta_0)$ is linear in $\delta_0$, and so can be represented as a matrix.
+The Lyapnuov exponents are defined from the $QR$ decomposition of $\varphi_{0,t}$: if
+\begin{equation}
+  \varphi_{0,t}=Q_tR_t
+\end{equation}
+where $Q_t$ is orthogonal, and $R_t$ is triangular with positive diagonal entries $(r_t^{(j)})_j$, then the Lyapunov exponents are
+\begin{equation}
+  \lim_{t\to\infty}\frac 1t\log|r_t^{(j)}|
+  .
+\end{equation}
+\bigskip
+
+\indent
+One problem to compute the Lyapunov exponents numerically is that the spectrum depends exponentially on time, and so has a tendency to blow up or shrink very rapidly, which leads to large truncation errors.
+To avoid these, we compute the tangent flow {\it in parts}: we split time into intervals:
+\begin{equation}
+  [0,t)=\bigcup_{i=0}^{N-1}[t_i,t_{i+1})
+  ,\quad
+  \varphi_{0,t}=\prod_{i=N-1}^0\varphi_{t_{i},t_{i+1}}
+  .
+\end{equation}
+At the thresholds between these intervals, we perform a QR decomposition:
+\begin{equation}
+  Q_0R_0:=\varphi_{0,t_0}
+  ,\quad
+  Q_iR_i:=\varphi_{t_{i-1},t_i}Q_{i-1}
+\end{equation}
+where $Q_i$ is orthogonal and $R_i$ is upper triangular with $\geqslant 0$ diagonal entries (in addition, we assume these are not $0$).
+Doing so, we find
+\begin{equation}
+  \varphi_{0,t}=Q_{N-1}\prod_{i=N-1}^0R_i
+\end{equation}
+and since the product of triangular matrices is triangular and the diagonal elements multiply, we find that the Lyapunov exponents are given by
+\begin{equation}
+  \lim_{t_{N-1}\to\infty}\frac1{t_{N-1}}\sum_{i=N-1}^0\log|r_i^{(j)}|
+  .
+\end{equation}
+To do the computation numerically, we drop the limit, and compute the logarithm of the product of the diagonal entries of $R_i$.
+\bigskip
+
+\indent
+In practice, we approximate $\varphi_{t_{i-1},t_i}$ by running a Runge-Kutta algorithm for the tangent flow equation.
+Because the tangent flow equation depends on $u(t)$, we must run it at times for which $u$ has been computed, so the Runge-Kutta algorithm for the tangent flow cannot be an adaptive step method, and should be one of {\tt RK4} for the fourth order algorithm or {\tt RK2} for the second order algorithm.
+To obtain the full matrix, we consider every element of the canonical basis as an initial condition $\delta_0$.
+We then iterate the Runge-Kutta algorithm until the time $t_0$ (chosen in one of two ways, see below), at which point we perform a QR decomposition, save the diagonal entries of $R$, replace the family of initial conditions with the columns of $Q$, and continue the flow from there.
+The choice of the times $t_i$ can be done either by fixed-length intervals, specified with the option {\tt lyapunov\_reset}, or the QR decomposition can be triggered whenever $\|\delta\|_1$ exceeds a threshold, specified in {\tt lyapunov\_maxu} (after all, the intervals are used to prevent $\delta$ from becoming too large). 
+\bigskip
+
+\indent
+To compute the Lyapunov exponents, we thus need the Jacobians of $\theta$ and $\psi$.
+Note that, by the linearity of the tangent flow equation,
+\begin{equation}
+  ((D \theta(\hat u))\delta)_{k}
+  =
+  \mathcal Re(Df(\hat u)(\delta_{\mathrm r}+i\delta_{\mathrm i}))
+  ,\quad
+  ((D \psi(\hat u))\delta)_{k}
+  =
+  \mathcal Im(Df(\hat u)(\delta_{\mathrm r}+i\delta_{\mathrm i}))
+  .
+\end{equation}
+For the irreversible equation,
+\begin{equation}
+  f(\hat u)=
+  -\frac{4\pi^2}{L^2}\nu k^2\hat u_k+\hat g_k
+  +\frac{4\pi^2}{L^2|k|}T(\hat u,k)
+\end{equation}
+and
+\begin{equation}
+  ((D f(\hat u))\delta)_{k}
+  =
+  -\frac{4\pi^2}{L^2}\nu k^2\delta_{k}
+  +\frac{4\pi^2}{L^2|k|}DT(\hat u,k)\delta
+\end{equation}
+\begin{equation}
+  DT(\hat u,k)\delta
+  =
+  \sum_{\displaystyle\mathop{\scriptstyle p,q\in\mathbb Z^2}_{p+q=k}}
+  \left(\frac{(q\cdot p^\perp)|q|}{|p|}+\frac{(p\cdot q^\perp)|p|}{|q|}\right)\hat u_p(\delta_{q,\mathrm r}+i\delta_{q,\mathrm i})
  .
 \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)
@@ -428,88 +575,97 @@ The parameter $\epsilon$ can be set using the parameter {\tt D\_epsilon}.
 %  \right)\hat u_{k-\ell}
 %  .
 %\end{equation}
-\bigskip
-
-\indent
-The Lyapunov exponents are then defined as
+For the reversible equation,
 \begin{equation}
-  \frac1{t_n}\log\mathrm{spec}(\pi_{t_n})
-  ,\quad
-  \pi_{t_n}:=\prod_{i=1}^nD\hat u(t_i)
-  .
-\end{equation}
-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.
-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:
-\begin{equation}
-  \pi_{\alpha\mathfrak L_r}
-  =R^{(\alpha)}Q^{(\alpha)}
-\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
-\begin{equation}
-  \frac1{\alpha\mathfrak L_r}\sum_{\beta=1}^\alpha\log\mathrm{spec}(Q^{(\beta)})
-  .
-\end{equation}
-
-\subsubsection{Numerical instability}.
-In order to prevent the algorithm from blowing up, it is necessary to impose the reality of $u(x)$ by hand, otherwise, truncation errors build up, and lead to divergences.
-It is sufficient to ensure that the convolution term $T(\hat u,k)$ satisfies $T(\hat u,-k)=T(\hat u,k)^*$.
-After imposing this condition, the algorithm no longer blows up, but it is still unstable (for instance, increasing $K_1$ or $K_2$ leads to very different results).
-
-\subsection{Reversible equation}
-\indent The reversible equation is similar to\-~(\ref{ins}) but instead of fixing the viscosity, we fix the enstrophy\-~\cite{Ga22}.
-It is defined directly in Fourier space:
-\begin{equation}
-  \partial_t\hat U_k=
-  -\frac{4\pi^2}{L^2}\alpha(\hat U) k^2\hat U_k+\hat G_k
-  -i\frac{2\pi}L\sum_{\displaystyle\mathop{\scriptstyle p,q\in\mathbb Z^2}_{p+q=k}}
-  (q\cdot\hat U_p)\hat U_q
-  ,\quad
-  k\cdot\hat U_k=0
-\end{equation}
-where $\alpha$ is chosen such that the enstrophy is constant.
-In terms of $\hat u$\-~(\ref{udef}), (\ref{gdef}), (\ref{T}):
-\begin{equation}
-  \partial_t\hat u_k=
+  f(\hat u)=
  -\frac{4\pi^2}{L^2}\alpha(\hat u) k^2\hat u_k
  +\hat g_k
  +\frac{4\pi^2}{L^2|k|}T(\hat u,k)
-  .
-  \label{rns_k}
 \end{equation}
-To compute $\alpha$, we use the constancy of the enstrophy:
+so
 \begin{equation}
-  \sum_{k\in\mathbb Z^2}k^2\hat U_k\cdot\partial_t\hat U_k
-  =0
-\end{equation}
-which, in terms of $\hat u$ is
-\begin{equation}
-  \sum_{k\in\mathbb Z^2}k^2\hat u_k^*\partial_t\hat u_k
-  =0
-\end{equation}
-that is
-\begin{equation}
-  \frac{4\pi^2}{L^2}\alpha(\hat u)\sum_{k\in\mathbb Z^2}k^4|\hat u_k|^2
+  ((D f(\hat u))\delta)_k
  =
-  \sum_{k\in\mathbb Z^2}k^2\hat u_k^*\hat g_k
-  +\frac{4\pi^2}{L^2}\sum_{k\in\mathbb Z^2}|k|\hat u_k^*T(\hat u,k)
+  -\frac{4\pi^2}{L^2}\alpha(\hat u) k^2\delta_k
+  -\frac{4\pi^2}{L^2}k^2\hat u_k D\alpha(\hat u)\delta
+  +\frac{4\pi^2}{L^2|k|}DT(\hat u,k)\delta
 \end{equation}
-and so
+where
 \begin{equation}
-  \alpha(\hat u)
-  =\frac{\frac{L^2}{4\pi^2}\sum_k k^2\hat u_k^*\hat g_k+\sum_k|k|\hat u_k^*T(\hat u,k)}{\sum_kk^4|\hat u_k|^2}
-  .
-  \label{alpha}
-\end{equation}
-Note that, by\-~(\ref{realu})-(\ref{realT}),
-\begin{equation}
-  \alpha(\hat u)\in\mathbb R
+  D\alpha(\hat u)\delta
+  =
+  \frac{\frac{L^2}{4\pi^2}\sum_k k^2\hat \delta_k^*\hat g_k+\sum_k|k|\hat (\delta_k^*T(\hat u,k)+\hat u_k^*DT(\hat u,k)\delta)}{\sum_kk^4|\hat u_k|^2}
+  -\alpha(\hat u)\frac{2\sum_kk^4\delta_k^*\hat u_k}{\sum_kk^4|\hat u_k|^2}
  .
 \end{equation}


+%We compute this Jacobian numerically using a finite difference, by computing
+%\begin{equation}
+%  (D\mathfrak F_{t_n})_{k,p}:=\frac1\epsilon
+%  \left(\begin{array}{cc}
+%    \phi_k(\hat u+\epsilon\delta_p)-\phi_k(\hat u)&\phi_k(\hat u+i\epsilon\delta_p)-\phi_k(\hat u)\\
+%    \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)
+%  .
+%\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
+%
+%\indent
+%The Lyapunov exponents are then defined as
+%\begin{equation}
+%  \frac1{t_n}\log\mathrm{spec}(\pi_{t_n})
+%  ,\quad
+%  \pi_{t_n}:=\prod_{i=1}^nD\hat u(t_i)
+%  .
+%\end{equation}
+%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.
+%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:
+%\begin{equation}
+%  \pi_{\alpha\mathfrak L_r}
+%  =R^{(\alpha)}Q^{(\alpha)}
+%\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
+%\begin{equation}
+%  \frac1{\alpha\mathfrak L_r}\sum_{\beta=1}^\alpha\log\mathrm{spec}(Q^{(\beta)})
+%  .
+%\end{equation}
+
+

 \vfill
 \eject
--- a/src/complex_tools.c
+++ b/src/complex_tools.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2017-2024 Ian Jauslin
+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.
--- a/src/complex_tools.h
+++ b/src/complex_tools.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2017-2024 Ian Jauslin
+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.
--- a/src/constants.cpp
+++ b/src/constants.cpp
@@ -1,5 +1,5 @@
 /*
-Copyright 2017-2024 Ian Jauslin
+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.
@@ -21,6 +21,7 @@ limitations under the License.
 #define COMMAND_QUIET 3
 #define COMMAND_RESUME 4
 #define COMMAND_LYAPUNOV 5
+#define COMMAND_ENSTROPHY_PRINT_INIT 6

 #define DRIVING_ZERO 1
 #define DRIVING_TEST 2
@@ -48,3 +49,6 @@ limitations under the License.

 #define FIX_ENSTROPHY 1
 #define FIX_ENERGY 2
+
+#define LYAPUNOV_TRIGGER_TIME 1
+#define LYAPUNOV_TRIGGER_SIZE 2
--- a/src/driving.c
+++ b/src/driving.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2017-2024 Ian Jauslin
+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.
--- a/src/driving.h
+++ b/src/driving.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2017-2024 Ian Jauslin
+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.
--- a/src/dstring.c
+++ b/src/dstring.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2017-2024 Ian Jauslin
+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.
--- a/src/dstring.h
+++ b/src/dstring.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2017-2024 Ian Jauslin
+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.
--- a/src/init.c
+++ b/src/init.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2017-2024 Ian Jauslin
+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.
--- a/src/init.h
+++ b/src/init.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2017-2024 Ian Jauslin
+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.
--- a/src/int_tools.c
+++ b/src/int_tools.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2017-2024 Ian Jauslin
+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.
--- a/src/int_tools.h
+++ b/src/int_tools.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2017-2024 Ian Jauslin
+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.
--- a/src/io.c
+++ b/src/io.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2017-2024 Ian Jauslin
+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.
@@ -51,6 +51,30 @@ int write_vec_bin(_Complex double* vec, int K1, int K2, FILE* file){
  return 0;
 }

+// write complex vector (stored as 2 doubles) indexed by k1,k2 to file in binary format
+int write_vec2_bin(double* vec, int K1, int K2, FILE* file){
+  // do nothing if there is no file
+  if(file==NULL){
+    return 0;
+  }
+
+  fwrite(vec, sizeof(double), 2*(K1*(2*K2+1)+K2), file);
+
+  return 0;
+}
+
+// write complex matrix (stored as 2 doubles) indexed by k1,k2 to file in binary format
+int write_mat2_bin(double* mat, int K1, int K2, FILE* file){
+  // do nothing if there is no file
+  if(file==NULL){
+    return 0;
+  }
+
+  fwrite(mat, sizeof(double), 4*(K1*(2*K2+1)+K2)*(K1*(2*K2+1)+K2), file);
+
+  return 0;
+}
+
 // read complex vector indexed by k1,k2 from file
 int read_vec(_Complex double* out, int K1, int K2, FILE* file){
  int kx,ky;
@@ -149,15 +173,53 @@ int read_vec(_Complex double* out, int K1, int K2, FILE* file){

 // read complex vector indexed by k1,k2 from file in binary format
 int read_vec_bin(_Complex double* out, int K1, int K2, FILE* file){
-  char c;
-  int ret;
-
  // do nothing if there is no file
  if(file==NULL){
    return 0;
  }

  // seek past initial comments
+  seek_past_initial_comments(file);
+
+  fread(out, sizeof(_Complex double), K1*(2*K2+1)+K2, file);
+
+  return 0;
+}
+
+// read complex vector (represented as 2 doubles) indexed by k1,k2 from file in binary format
+int read_vec2_bin(double* out, int K1, int K2, FILE* file){
+  if(file==NULL){
+    return 0;
+  }
+
+  // seek past initial comments
+  seek_past_initial_comments(file);
+
+  fread(out, sizeof(double), 2*(K1*(2*K2+1)+K2), file);
+  return 0;
+}
+
+// read complex matrix (represented as 2 doubles) indexed by k1,k2 from file in binary format
+int read_mat2_bin(double* out, int K1, int K2, FILE* file){
+  if(file==NULL){
+    return 0;
+  }
+
+  // seek past initial comments
+  seek_past_initial_comments(file);
+
+  fread(out, sizeof(double), 4*(K1*(2*K2+1)+K2)*(K1*(2*K2+1)+K2), file);
+  return 0;
+}
+
+// ignore comments at beginning of file
+int seek_past_initial_comments(FILE* file){
+  char c;
+  int ret;
+
+  if(file==NULL){
+    return 0;
+  }
  while(true){
    ret=fscanf(file, "%c", &c);
    if (ret==1 && c=='#'){
@@ -184,8 +246,6 @@ int read_vec_bin(_Complex double* out, int K1, int K2, FILE* file){
    }
  }

-  fread(out, sizeof(_Complex double), K1*(2*K2+1)+K2, file);
-
  return 0;
 }

@@ -198,23 +258,38 @@ int remove_entry(
  char* rw_ptr;
  char* bfr;
  char* entry_ptr=entry;
+  // whether to write the entry
  int go=1;
+  // whether the pointer is at the beginning of the entry
+  int at_top=1;

  for(ptr=param_str, rw_ptr=ptr; *ptr!='\0'; ptr++){
-    for(bfr=ptr,entry_ptr=entry; *bfr==*entry_ptr; bfr++, entry_ptr++){
-      // check if reached end of entry
-      if(*(bfr+1)=='=' && *(entry_ptr+1)=='\0'){
-	go=0;
-	break;
+    // only match entries if one is at the beginning of an entry
+    if(at_top==1){
+      // check that the entry under ptr matches entry
+      for(bfr=ptr,entry_ptr=entry; *bfr==*entry_ptr; bfr++, entry_ptr++){
+	// check if reached end of entry
+	if(*(bfr+1)=='=' && *(entry_ptr+1)=='\0'){
+	  // match: do not write entry
+	  go=0;
+	  break;
+	}
      }
    }
+
+    // write entry
    if(go==1){
      *rw_ptr=*ptr;
      rw_ptr++;
    }
+
+    // next iterate will no longer be at the beginning of the entry
+    at_top=0;
+
    //reset
    if(*ptr==';'){
      go=1;
+      at_top=1;
    }
  }
  *rw_ptr='\0';
@@ -262,6 +337,8 @@ int save_state(
    strcpy(params, params_string);
    remove_entry(params, "starting_time");
    remove_entry(params, "init");
+    remove_entry(params, "init_enstrophy");
+    remove_entry(params, "init_energy");
    if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){
      remove_entry(params, "delta");
    }
--- a/src/io.h
+++ b/src/io.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2017-2024 Ian Jauslin
+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.
@@ -22,10 +22,21 @@ limitations under the License.
 // write complex vector indexed by k1,k2 to file
 int write_vec(_Complex double* u, int K1, int K2, FILE* file);
 int write_vec_bin(_Complex double* u, int K1, int K2, FILE* file);
+// write complex vector (stored as 2 doubles) indexed by k1,k2 to file in binary format
+int write_vec2_bin(double* vec, int K1, int K2, FILE* file);
+// write complex matrix (stored as 2 doubles) indexed by k1,k2 to file in binary format
+int write_mat2_bin(double* mat, int K1, int K2, FILE* file);

 // read complex vector indexed by k1,k2 from file
 int read_vec(_Complex double* u, int K1, int K2, FILE* file);
 int read_vec_bin(_Complex double* u, int K1, int K2, FILE* file);
+// read complex vector (represented as 2 doubles) indexed by k1,k2 from file in binary format
+int read_vec2_bin(double* out, int K1, int K2, FILE* file);
+// read complex matrix (represented as 2 doubles) indexed by k1,k2 from file in binary format
+int read_mat2_bin(double* out, int K1, int K2, FILE* file);
+
+// ignore comments at beginning of file
+int seek_past_initial_comments(FILE* file);

 // remove an entry from params string (inplace)
 int remove_entry(char* param_str, char* entry);
--- a/src/lyapunov.c
+++ b/src/lyapunov.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2017-2024 Ian Jauslin
+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.
@@ -16,6 +16,7 @@ limitations under the License.

 #include "constants.cpp"
 #include "lyapunov.h"
+#include "io.h"
 #include <openblas64/cblas.h>
 #include <openblas64/lapacke.h>
 #include <math.h>
@@ -32,120 +33,130 @@ 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
+  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* 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<MATSIZE;i++){
-    for(j=0;j<MATSIZE;j++){
-      Du_prod[i*MATSIZE+j]=0.;
+  // initialize flow and averages
+  for (i=0;i<MATSIZE;i++){
+    for (j=0;j<MATSIZE;j++){
+      flow[i*MATSIZE+j]=flow0[i*MATSIZE+j];
    }
-  }
-  for(i=0;i<MATSIZE;i++){
-    Du_prod[i*MATSIZE+i]=1.;
+    lyapunov_avg[i]=lyapunov_avg0[i];
  }

-  // store step (useful for adaptive step methods
-  double prev_step=delta;
+  // store step (useful for adaptive step methods)
  double step=delta;
  double next_step=step;

-  // init average
-  for (i=0; i<MATSIZE; i++){
-    lyapunov_avg[i]=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_dgemm(CblasColMajor, CblasNoTrans, CblasTrans, MATSIZE, MATSIZE, MATSIZE, 1., Du_prod, MATSIZE, Du, MATSIZE, 0., tmp10, MATSIZE);
-    // copy back to Du_prod
-    double* move=tmp10;
-    tmp10=Du_prod;
-    Du_prod=move;
+    // compute u first
+    // if using an adaptive step, the step for the tangent flow is set by this computation of u
+    // use fftu1 as a tmp fft vector (it hasn't been used yet this step)
+    ns_step(algorithm, u, K1, K2, N1, N2, nu, &step, &next_step, adaptive_tolerance, adaptive_factor, max_delta, adaptive_cost, L, g, time, starting_time, fft1, fftu1, ifft, &tmp4, &tmp5, tmp6, tmp7, tmp8, tmp9, tmp10, irreversible, keep_en_cst, target_en);
+    // compute tangent flow
+    lyapunov_D_step(flow, u, K1, K2, N1, N2, nu, step, algorithm_lyapunov, L, g, fftu1, fftu2, fftu3, fftu4, fft1, ifft, tmp1, tmp2, tmp3, tmp4, irreversible);

    // increment time
    time+=step;

-    // reset Jacobian
-    if(time>(n+1)*lyapunov_reset){
+    double norm=0;
+    if(lyapunov_trigger==LYAPUNOV_TRIGGER_SIZE){
+      // size of flow (for reset)
+      for(i=0;i<MATSIZE;i++){
+	for(j=0;j<MATSIZE;j++){
+	  norm+=flow[i*MATSIZE+j]*flow[i*MATSIZE+j]/MATSIZE;
+	  if(sqrt(norm)>lyapunov_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_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; i<MATSIZE; i++){
-	lyapunov[i]=log(fabs(Du_prod[i*MATSIZE+i]))/(time-prevtime);
+	lyapunov[i]=log(fabs(flow[i*MATSIZE+i]))/(time-prevtime);
      }

-      // sort lyapunov exponents
-      qsort(lyapunov, MATSIZE, sizeof(double), compare_double);
+      //// sort lyapunov exponents
+      //qsort(lyapunov, MATSIZE, sizeof(double), compare_double);

      // average lyapunov
      for(i=0; i<MATSIZE; i++){
 	// exclude inf
-	if((! isinf(lyapunov[i])) && (time>starting_time)){
-	  lyapunov_avg[i]=lyapunov_avg[i]*(prevtime-starting_time)/(time-starting_time)+lyapunov[i]*(time-prevtime)/(time-starting_time);
+	if((! isinf(lyapunov[i])) && (time>lyapunov_startingtime)){
+	  lyapunov_avg[i]=lyapunov_avg[i]*(prevtime-lyapunov_startingtime)/(time-lyapunov_startingtime)+lyapunov[i]*(time-prevtime)/(time-lyapunov_startingtime);
 	}
      }

@@ -154,24 +165,43 @@ int lyapunov(
 	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]);
-      }
+      fprintf(stderr,"% .15e\n",time);
+      //// print largest and smallest lyapunov exponent to stderr
+      //if(MATSIZE>0){
+      //  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;
    }
+
+    // catch abort signal
+    if (g_abort){
+      // print u to stderr if no savefile
+      if (savefile==NULL){
+	savefile=stderr;
+      }
+      break;
+    }
  }

-  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);
+  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) {
@@ -183,168 +213,603 @@ 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;
+  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++){
-      // derivative in the real direction
-      // finite difference vector
-      for (i=0;i<USIZE;i++){
-	if(i==klookup_sym(lx,ly,K2)){
-	  tmp1[i]=un[i]+epsilon;
-	}else{
-	  tmp1[i]=un[i];
-	}
-      }
-      // compute step 
-      // reinitialize step
-      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);
-      // derivative
-      for (i=0;i<USIZE;i++){
-	// use row major order
-	out[2*klookup_sym(lx,ly,K2)*MATSIZE+2*i]=(__real__ (tmp1[i]-un_next[i]))/epsilon;
-	out[2*klookup_sym(lx,ly,K2)*MATSIZE+2*i+1]=(__imag__ (tmp1[i]-un_next[i]))/epsilon;
-      }
+      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);
+  }

-      // derivative in the imaginary direction
-      // finite difference vector
-      for (i=0;i<USIZE;i++){
-	if(i==klookup_sym(lx,ly,K2)){
-	  tmp1[i]=un[i]+epsilon*I;
-	}else{
-	  tmp1[i]=un[i];
+  // 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);
 	}
      }
-      // compute step 
-      // reinitialize step
-      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);
-      // compute derivative
-      for (i=0;i<USIZE;i++){
-	// use row major order
-	out[(2*klookup_sym(lx,ly,K2)+1)*MATSIZE+2*i]=(__real__ (tmp1[i]-un_next[i]))/epsilon;
-	out[(2*klookup_sym(lx,ly,K2)+1)*MATSIZE+2*i+1]=(__imag__ (tmp1[i]-un_next[i]))/epsilon;
-      }
    }
  }

  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;i<MATSIZE;i++){
+    tmp3[i]=flow[i]+delta/6*tmp1[i];
+  }
+
+  // d+h*k1/2
+  for(i=0;i<MATSIZE;i++){
+    tmp2[i]=flow[i]+delta/2*tmp1[i];
+  }
+  // k2
+  lyapunov_D_rhs(tmp1, tmp2, u, K1, K2, N1, N2, nu, L, g, T, fftu1, fftu2, fftu3, fftu4, fft1, ifft, irreversible);
+  // add to output
+  for(i=0;i<MATSIZE;i++){
+    tmp3[i]+=delta/3*tmp1[i];
+  }
+
+  // d+h*k2/2
+  for(i=0;i<MATSIZE;i++){
+    tmp2[i]=flow[i]+delta/2*tmp1[i];
+  }
+  // k3
+  lyapunov_D_rhs(tmp1, tmp2, u, K1, K2, N1, N2, nu, L, g, T, fftu1, fftu2, fftu3, fftu4, fft1, ifft, irreversible);
+  // add to output
+  for(i=0;i<MATSIZE;i++){
+    tmp3[i]+=delta/3*tmp1[i];
+  }
+
+  // d+h*k3
+  for(i=0;i<MATSIZE;i++){
+    tmp2[i]=flow[i]+delta*tmp1[i];
+  }
+  // k4
+  lyapunov_D_rhs(tmp1, tmp2, u, K1, K2, N1, N2, nu, L, g, T, fftu1, fftu2, fftu3, fftu4, fft1, ifft, irreversible);
+  // add to output
+  for(i=0;i<MATSIZE;i++){
+    flow[i]=tmp3[i]+delta/6*tmp1[i];
+  }
+
+  return(0);
+}
+
+// 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
+){
+  int i;
+
+  // k1
+  lyapunov_D_rhs(tmp1, flow, u, K1, K2, N1, N2, nu, L, g, T, fftu1, fftu2, fftu3, fftu4, fft1, ifft, irreversible);
+
+  // u+h*k1/2
+  for(i=0;i<MATSIZE;i++){
+    tmp2[i]=flow[i]+delta/2*tmp1[i];
+  }
+  // k2
+  lyapunov_D_rhs(tmp1, tmp2, u, K1, K2, N1, N2, nu, L, g, T, fftu1, fftu2, fftu3, fftu4, fft1, ifft, irreversible);
+  // add to output
+  for(i=0;i<MATSIZE;i++){
+    flow[i]+=delta*tmp1[i];
+  }
+
+  return(0);
+}
+
+// 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
+){
+  int kx,ky;
+  int i;
+
+  // compute DT
+  lyapunov_D_T(flow,K1,K2,N1,N2,fftu1,fftu2,fftu3,fftu4,fft1,ifft);
+
+  for(i=0; i<K1*(2*K2+1)+K2; i++){
+    out[i]=0;
+  }
+  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*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;i<N1*N2;i++){
+    ifft.fft[i]=fftu1.fft[i]*fftu3.fft[i]-fftu2.fft[i]*fftu4.fft[i];
+  }
+  // inverse fft
+  fftw_execute(ifft.fft_plan);
+
+  return(0);
+}
+
+// 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
+){
+  int i;
+  int kx,ky;
+
+  // F(px/|p|*u)*F(qy*|q|*delta)
+  // init to 0
+  for(i=0; i<N1*N2; i++){
+    fft1.fft[i]=0;
+  }
+  // fill modes
+  for(kx=-K1;kx<=K1;kx++){
+    for(ky=-K2;ky<=K2;ky++){
+      if(kx!=0 || ky!=0){
+        fft1.fft[klookup(kx,ky,N1,N2)]=ky*sqrt(kx*kx+ky*ky)*lyapunov_delta_getval_sym(flow, kx,ky,K2)/N2;
+      }
+    }
+  }
+  // fft
+  fftw_execute(fft1.fft_plan);
+  // write to ifft
+  for(i=0;i<N1*N2;i++){
+    ifft.fft[i]=fftu1.fft[i]*fft1.fft[i];
+  }
+
+  // F(px/|p|*delta)*F(qy*|q|*p)
+  // init to 0
+  for(i=0; i<N1*N2; i++){
+    fft1.fft[i]=0;
+  }
+  // fill modes
+  for(kx=-K1;kx<=K1;kx++){
+    for(ky=-K2;ky<=K2;ky++){
+      if(kx!=0 || ky!=0){
+        fft1.fft[klookup(kx,ky,N1,N2)]=kx/sqrt(kx*kx+ky*ky)*lyapunov_delta_getval_sym(flow, kx,ky,K2)/N1;
+      }
+    }
+  }
+  // fft
+  fftw_execute(fft1.fft_plan);
+  // write to ifft
+  for(i=0;i<N1*N2;i++){
+    ifft.fft[i]=fftu2.fft[i]*fft1.fft[i];
+  }
+
+  // F(py/|p|*u)*F(qx*|q|*delta)
+  // init to 0
+  for(i=0; i<N1*N2; i++){
+    fft1.fft[i]=0;
+  }
+  // fill modes
+  for(kx=-K1;kx<=K1;kx++){
+    for(ky=-K2;ky<=K2;ky++){
+      if(kx!=0 || ky!=0){
+        fft1.fft[klookup(kx,ky,N1,N2)]=kx*sqrt(kx*kx+ky*ky)*lyapunov_delta_getval_sym(flow, kx,ky,K2)/N2;
+      }
+    }
+  }
+  // fft
+  fftw_execute(fft1.fft_plan);
+  // write to ifft
+  for(i=0;i<N1*N2;i++){
+    ifft.fft[i]=ifft.fft[i]-fftu3.fft[i]*fft1.fft[i];
+  }
+
+  // F(py/|p|*delta)*F(qx*|q|*u)
+  // init to 0
+  for(i=0; i<N1*N2; i++){
+    fft1.fft[i]=0;
+  }
+  // fill modes
+  for(kx=-K1;kx<=K1;kx++){
+    for(ky=-K2;ky<=K2;ky++){
+      if(kx!=0 || ky!=0){
+        fft1.fft[klookup(kx,ky,N1,N2)]=ky/sqrt(kx*kx+ky*ky)*lyapunov_delta_getval_sym(flow, kx,ky,K2)/N1;
+      }
+    }
+  }
+  // fft
+  fftw_execute(fft1.fft_plan);
+  // write to ifft
+  for(i=0;i<N1*N2;i++){
+    ifft.fft[i]=ifft.fft[i]-fftu4.fft[i]*fft1.fft[i];
+  }
+
+  // inverse fft
+  fftw_execute(ifft.fft_plan);
+  
+  return 0;
+}
+
+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
+){
+  int kx,ky;
+
+  _Complex double num=0.;
+  _Complex double denom=0.;
+
+  for(kx=0;kx<=K1;kx++){
+    for(ky=(kx>0 ? -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; i<N1*N2; i++){
+    fftu1.fft[i]=0;
+    fftu2.fft[i]=0;
+    fftu3.fft[i]=0;
+    fftu4.fft[i]=0;
+  }
+  // fill modes
+  for(kx=-K1;kx<=K1;kx++){
+    for(ky=-K2;ky<=K2;ky++){
+      if(kx!=0 || ky!=0){
+        fftu1.fft[klookup(kx,ky,N1,N2)]=kx/sqrt(kx*kx+ky*ky)*getval_sym(u, kx,ky,K2)/N1;
+        fftu2.fft[klookup(kx,ky,N1,N2)]=ky*sqrt(kx*kx+ky*ky)*getval_sym(u, kx,ky,K2)/N2;
+        fftu3.fft[klookup(kx,ky,N1,N2)]=ky/sqrt(kx*kx+ky*ky)*getval_sym(u, kx,ky,K2)/N1;
+        fftu4.fft[klookup(kx,ky,N1,N2)]=kx*sqrt(kx*kx+ky*ky)*getval_sym(u, kx,ky,K2)/N2;
+      }
+    }
+  }
+
+  // fft
+  fftw_execute(fftu1.fft_plan);
+  fftw_execute(fftu2.fft_plan);
+  fftw_execute(fftu3.fft_plan);
+  fftw_execute(fftu4.fft_plan);
+
+  return(0);
+}

 int lyapunov_init_tmps(
  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 ** 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,
+  _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,
  int K1,
  int K2,
  int N1,
  int N2,
  unsigned int nthreads,
-  unsigned int algorithm
+  unsigned int algorithm,
+  unsigned int algorithm_lyapunov
 ){
-  ns_init_tmps(u, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9, fft1, fft2, ifft, K1, K2, N1, N2, nthreads, algorithm);
+  // allocate tmp vectors for computation
+  if(algorithm_lyapunov==ALGORITHM_RK2){
+    *tmp1=calloc(MATSIZE,sizeof(double));
+    *tmp2=calloc(MATSIZE,sizeof(double));
+  } else if (algorithm_lyapunov==ALGORITHM_RK4){
+    *tmp1=calloc(MATSIZE,sizeof(double));
+    *tmp2=calloc(MATSIZE,sizeof(double));
+    *tmp3=calloc(MATSIZE,sizeof(double));
+  } else {
+    fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm_lyapunov);
+  };
+  *tmp11=calloc(MATSIZE*MATSIZE,sizeof(double));
+
+  ns_init_tmps(u, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9, tmp10, fft1, fftu1, ifft, K1, K2, N1, N2, nthreads, algorithm);
+
+  // prepare vectors for fft
+  fftu2->fft=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);
 }
+
+// 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;
+}
--- a/src/lyapunov.h
+++ b/src/lyapunov.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2017-2024 Ian Jauslin
+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.
@@ -20,18 +20,44 @@ 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_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, double lyapunov_startingtime, FILE* savefile, FILE* utfile, const char* cmd_string, const char* params_string, const char* savefile_string, const char* utfile_string);

 // 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);
+
+// 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);

 #endif

--- a/src/main.c
+++ b/src/main.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2017-2024 Ian Jauslin
+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.
@@ -28,6 +28,7 @@ limitations under the License.
 #include "dstring.h"
 #include "init.h"
 #include "int_tools.h"
+#include "io.h"
 #include "lyapunov.h"
 #include "navier-stokes.h"

@@ -56,11 +57,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;
-  double D_epsilon;
+  unsigned int lyapunov_trigger;
  bool print_alpha;
 } nstrophy_parameters;

@@ -81,6 +83,8 @@ int args_from_file(dstring* params, unsigned int* command, unsigned int* nthread
 _Complex double* set_driving(nstrophy_parameters parameters);
 // set initial condition
 _Complex double* set_init(nstrophy_parameters* parameters);
+// set initial tangent flow for lyapunov exponents
+int set_lyapunov_flow_init( double** lyapunov_flow0, double** lyapunov_avg0, double* lyapunov_prevtime, double* lyapunov_startingtime, bool fromfile, nstrophy_parameters parameters);

 // signal handler
 void sig_handler (int signo);
@@ -114,6 +118,10 @@ int main (
  unsigned int nthreads=1;
  _Complex double* u0;
  _Complex double *g;
+  double* lyapunov_flow0;
+  double* lyapunov_avg0;
+  double lyapunov_prevtime;
+  double lyapunov_startingtime;
  dstring savefile_str;
  dstring utfile_str;
  dstring initfile_str;
@@ -121,6 +129,7 @@ int main (
  dstring resumefile_str;
  FILE* savefile=NULL;
  FILE* utfile=NULL;
+  bool resume=false;

  command=0;

@@ -160,6 +169,8 @@ int main (

  // if command is 'resume', then read args from file
  if(command==COMMAND_RESUME){
+    // remember that the original command was resume (to set values from init file)
+    resume=true;
    ret=args_from_file(&param_str, &command, &nthreads, &savefile_str, &utfile_str, dstring_to_str_noinit(&resumefile_str));
    if(ret<0){
      dstring_free(param_str);
@@ -223,6 +234,20 @@ int main (
  g=set_driving(parameters);
  // set initial condition
  u0=set_init(&parameters);
+  // read extra values from init file when resuming a computation
+  if(resume){
+    // read start time
+    fread(&(parameters.starting_time), sizeof(double), 1, parameters.initfile);
+    // if adaptive step algorithm
+    if(parameters.algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){
+      // read delta
+      fread(&(parameters.delta), sizeof(double), 1, parameters.initfile);
+    }
+  }
+  // set initial condition for the lyapunov flow
+  if (command==COMMAND_LYAPUNOV){
+    set_lyapunov_flow_init(&lyapunov_flow0, &lyapunov_avg0, &lyapunov_prevtime, &lyapunov_startingtime, resume, parameters);
+  }

  // if init_enstrophy is not set in the parameters, then compute it from the initial condition
  if (parameters.init_enstrophy_or_energy!=FIX_ENSTROPHY){
@@ -254,6 +279,10 @@ int main (
      dstring_free(drivingfile_str);
      free(g);
      free(u0);
+      if (command==COMMAND_LYAPUNOV){
+	free(lyapunov_flow0);
+	free(lyapunov_avg0);
+      }
      return(-1);
    }
  }
@@ -270,6 +299,10 @@ int main (
      dstring_free(drivingfile_str);
      free(g);
      free(u0);
+      if (command==COMMAND_LYAPUNOV){
+	free(lyapunov_flow0);
+	free(lyapunov_avg0);
+      }
      return(-1);
    }
  }
@@ -284,7 +317,7 @@ int main (

  // run command
  if (command==COMMAND_UK){
-    uk(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, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_enstrophy, parameters.algorithm, parameters.print_freq, parameters.starting_time, nthreads, savefile);
+    uk(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, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_enstrophy, parameters.algorithm, parameters.print_freq, parameters.starting_time, nthreads, savefile, utfile, (char*)argv[0], dstring_to_str_noinit(&param_str), dstring_to_str_noinit(&savefile_str), dstring_to_str_noinit(&utfile_str));
  }
  else if(command==COMMAND_ENSTROPHY){
    // register signal handler to handle aborts
@@ -292,11 +325,17 @@ int main (
    signal(SIGTERM, sig_handler);
    enstrophy(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, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_enstrophy, parameters.algorithm, parameters.print_freq, parameters.starting_time, parameters.print_alpha, nthreads, savefile, utfile, (char*)argv[0], dstring_to_str_noinit(&param_str), dstring_to_str_noinit(&savefile_str), dstring_to_str_noinit(&utfile_str));
  }
+  else if(command==COMMAND_ENSTROPHY_PRINT_INIT){
+    enstrophy_print_init(parameters.K1, parameters.K2, parameters.L, u0, g);
+  }
  else if(command==COMMAND_QUIET){
    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);
+    // register signal handler to handle aborts
+    signal(SIGINT, sig_handler);
+    signal(SIGTERM, sig_handler);
+    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, lyapunov_flow0, lyapunov_avg0, lyapunov_prevtime, lyapunov_startingtime, savefile, utfile, (char*)argv[0], dstring_to_str_noinit(&param_str), dstring_to_str_noinit(&savefile_str), dstring_to_str_noinit(&utfile_str));
  }
  else if(command==0){
    fprintf(stderr, "error: no command specified\n");
@@ -305,6 +344,10 @@ int main (

  free(g);
  free(u0);
+  if (command==COMMAND_LYAPUNOV){
+    free(lyapunov_flow0);
+    free(lyapunov_avg0);
+  }

  // free strings
  dstring_free(param_str);
@@ -521,6 +564,9 @@ int read_args(
      else if(strcmp(argv[i], "enstrophy")==0){
 	*command=COMMAND_ENSTROPHY;
      }
+      else if(strcmp(argv[i], "enstrophy_print_init")==0){
+	*command=COMMAND_ENSTROPHY_PRINT_INIT;
+      }
      else if(strcmp(argv[i], "quiet")==0){
 	*command=COMMAND_QUIET;
      }
@@ -539,7 +585,7 @@ int read_args(
  }

  // check that if the command is 'resume', then resumefile has been set
-  if(*command==COMMAND_RESUME && resumefile_str->length==0){
+  if(*command==COMMAND_RESUME && (resumefile_str==NULL || resumefile_str->length==0)){
    fprintf(stderr, "error: 'resume' command used, but no resume file\n");
    print_usage();
    return(-1);
@@ -576,6 +622,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 +697,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,20 +896,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){
-      fprintf(stderr, "error: parameter 'D_epsilon' should be a double\n       got '%s'\n",rhs);
-      return(-1);
-    }
-  }
  else if (strcmp(lhs,"driving")==0){
    if (strcmp(rhs,"zero")==0){
      parameters->driving=DRIVING_ZERO;
@@ -940,6 +981,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);
@@ -1055,12 +1137,6 @@ _Complex double* set_init(

  case INIT_FILE:
    init_file_bin(u0, parameters->K1, parameters->K2, parameters->initfile);
-    // read start time
-    fread(&(parameters->starting_time), sizeof(double), 1, parameters->initfile);
-    if(parameters->algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){
-      // read delta
-      fread(&(parameters->delta), sizeof(double), 1, parameters->initfile);
-    }
    break;

  case INIT_FILE_TXT:
@@ -1095,3 +1171,49 @@ _Complex double* set_init(
  
  return u0;
 }
+
+// set initial tangent flow for lyapunov exponents
+int set_lyapunov_flow_init(
+  double** lyapunov_flow0,
+  double** lyapunov_avg0,
+  double* lyapunov_prevtime,
+  double* lyapunov_startingtime,
+  bool fromfile, // whether to initialize from file
+  nstrophy_parameters parameters
+){
+  *lyapunov_flow0=calloc(4*(parameters.K1*(2*parameters.K2+1)+parameters.K2)*(parameters.K1*(2*parameters.K2+1)+parameters.K2),sizeof(double));
+  *lyapunov_avg0=calloc(2*(parameters.K1*(2*parameters.K2+1)+parameters.K2),sizeof(double));
+
+  // read from file or init from identity matrix
+  if(fromfile){
+    // read flow
+    read_mat2_bin(*lyapunov_flow0, parameters.K1, parameters.K2, parameters.initfile);
+    // read time of last QR decomposition
+    fread(lyapunov_prevtime, sizeof(double), 1, parameters.initfile);
+    // read time at which lyapunov computation was started
+    fread(lyapunov_startingtime, sizeof(double), 1, parameters.initfile);
+    // read averages
+    read_vec2_bin(*lyapunov_avg0, parameters.K1, parameters.K2, parameters.initfile);
+  } else {
+    // init with identity
+    int i,j;
+    for (i=0;i<2*(parameters.K1*(2*parameters.K2+1)+parameters.K2);i++){
+      for (j=0;j<2*(parameters.K1*(2*parameters.K2+1)+parameters.K2);j++){
+	if(i!=j){
+	  (*lyapunov_flow0)[i*2*(parameters.K1*(2*parameters.K2+1)+parameters.K2)+j]=0.;
+	} else {
+	  (*lyapunov_flow0)[i*2*(parameters.K1*(2*parameters.K2+1)+parameters.K2)+j]=1.;
+	}
+      }
+    }
+    // init prevtime
+    *lyapunov_prevtime=parameters.starting_time;
+    // init starting_time
+    *lyapunov_startingtime=parameters.starting_time;
+    // init averages
+    for(i=0;i<2*(parameters.K1*(2*parameters.K2+1)+parameters.K2);i++){
+      (*lyapunov_avg0)[i]=0;
+    }
+  }
+  return 0;
+}
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -1,5 +1,5 @@
 /*
-Copyright 2017-2024 Ian Jauslin
+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.
@@ -23,6 +23,8 @@ limitations under the License.
 #include <stdlib.h>
 #include <string.h>

+#define USIZE (K1*(2*K2+1)+K2)
+
 // compute solution as a function of time
 int uk(
  int K1,
@@ -46,7 +48,13 @@ int uk(
  double print_freq,
  double starting_time,
  unsigned int nthreads,
-  FILE* savefile
+  FILE* savefile,
+  FILE* utfile,
+  // for interrupt recovery
+  const char* cmd_string,
+  const char* params_string,
+  const char* savefile_string,
+  const char* utfile_string
 ){
  _Complex double* u;
  _Complex double* tmp1;
@@ -82,7 +90,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;

@@ -93,7 +101,6 @@ int uk(
    ns_step(algorithm, u, K1, K2, N1, N2, nu, &step, &next_step, adaptive_tolerance, adaptive_factor, max_delta, adaptive_cost, L, g, time, starting_time, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en);

    time+=step;
-    step=next_step;
    
    if(time>(n+1)*print_freq){
      n++;
@@ -113,11 +120,25 @@ int uk(
      fprintf(stderr,"\n");
      printf("\n");
    }
+
+    // get ready for next step
+    step=next_step;
+
+    // catch abort signal
+    if (g_abort){
+      break;
+    }
  }

  // TODO: update handling of savefile
-  // save final entry to savefile
-  write_vec_bin(u, K1, K2, savefile);
+  if(savefile!=NULL){
+    save_state(u, savefile, K1, K2, cmd_string, params_string, savefile_string, utfile_string, utfile, COMMAND_UK, algorithm, step, time, nthreads);
+  }
+
+  // save final u to utfile in txt format
+  if(utfile!=NULL){
+    write_vec(u, K1, K2, utfile);
+  }

  ns_free_tmps(u, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, fft1, fft2, ifft, algorithm);
  return(0);
@@ -223,10 +244,18 @@ int enstrophy(
      // print to stderr so user can follow along
      if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){
 	fprintf(stderr,"% .8e % .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",time, avg_a, avg_energy, avg_enstrophy, alpha, energy, enstrophy, step);
-	printf("% .15e % .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",time, avg_a, avg_energy, avg_enstrophy, alpha, energy, enstrophy, step);
+	if(K1>=1 && K2>=2){
+	  printf("% .15e % .15e % .15e % .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",time, avg_a, avg_energy, avg_enstrophy, alpha, energy, enstrophy, step, __real__ u[klookup_sym(1,1,K2)], __real__ u[klookup_sym(1,2,K2)]);
+	}else{
+	  printf("% .15e % .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",time, avg_a, avg_energy, avg_enstrophy, alpha, energy, enstrophy, step);
+	}
      } else {
 	fprintf(stderr,"% .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",time, avg_a, avg_energy, avg_enstrophy, alpha, energy, enstrophy);
-	printf("% .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",time, avg_a, avg_energy, avg_enstrophy, alpha, energy, enstrophy);
+	if(K1>=1 && K2>=2){
+	  printf("% .15e % .15e % .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",time, avg_a, avg_energy, avg_enstrophy, alpha, energy, enstrophy, __real__ u[klookup_sym(1,1,K2)], __real__ u[klookup_sym(1,2,K2)]);
+	}else{
+	  printf("% .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",time, avg_a, avg_energy, avg_enstrophy, alpha, energy, enstrophy);
+	}
      }

      // reset averages
@@ -268,6 +297,25 @@ int enstrophy(
  return(0);
 }

+// compute enstrophy, alpha for the initial condition (useful for debugging)
+int enstrophy_print_init(
+  int K1,
+  int K2,
+  double L,
+  _Complex double* u0,
+  _Complex double* g
+){
+  double alpha, enstrophy, energy;
+
+  alpha=compute_alpha(u0, K1, K2, g, L);
+  enstrophy=compute_enstrophy(u0, K1, K2, L);
+  energy=compute_energy(u0, K1, K2);
+
+  // print
+  printf("% .15e % .15e % .15e\n", alpha, energy, enstrophy);
+  return(0);
+}
+
 // compute solution as a function of time, but do not print anything (useful for debugging)
 int quiet(
  int K1,
@@ -353,30 +401,30 @@ int ns_init_tmps(
  unsigned int algorithm
 ){
  // velocity field
-  *u=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
+  *u=calloc(USIZE,sizeof(_Complex double));

  // allocate tmp vectors for computation
  if(algorithm==ALGORITHM_RK2){
-    *tmp1=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
-    *tmp2=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
+    *tmp1=calloc(USIZE,sizeof(_Complex double));
+    *tmp2=calloc(USIZE,sizeof(_Complex double));
  } else if (algorithm==ALGORITHM_RK4){
-    *tmp1=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
-    *tmp2=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
-    *tmp3=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
+    *tmp1=calloc(USIZE,sizeof(_Complex double));
+    *tmp2=calloc(USIZE,sizeof(_Complex double));
+    *tmp3=calloc(USIZE,sizeof(_Complex double));
  } else if (algorithm==ALGORITHM_RKF45 || algorithm==ALGORITHM_RKDP54){
-    *tmp1=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
-    *tmp2=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
-    *tmp3=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
-    *tmp4=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
-    *tmp5=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
-    *tmp6=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
-    *tmp7=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
+    *tmp1=calloc(USIZE,sizeof(_Complex double));
+    *tmp2=calloc(USIZE,sizeof(_Complex double));
+    *tmp3=calloc(USIZE,sizeof(_Complex double));
+    *tmp4=calloc(USIZE,sizeof(_Complex double));
+    *tmp5=calloc(USIZE,sizeof(_Complex double));
+    *tmp6=calloc(USIZE,sizeof(_Complex double));
+    *tmp7=calloc(USIZE,sizeof(_Complex double));
  } else if (algorithm==ALGORITHM_RKBS32){
-    *tmp1=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
-    *tmp2=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
-    *tmp3=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
-    *tmp4=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
-    *tmp5=calloc(K1*(2*K2+1)+K2,sizeof(_Complex double));
+    *tmp1=calloc(USIZE,sizeof(_Complex double));
+    *tmp2=calloc(USIZE,sizeof(_Complex double));
+    *tmp3=calloc(USIZE,sizeof(_Complex double));
+    *tmp4=calloc(USIZE,sizeof(_Complex double));
+    *tmp5=calloc(USIZE,sizeof(_Complex double));
  } else {
    fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
  };
@@ -462,7 +510,7 @@ int copy_u(
 ){
  int i;

-  for(i=0;i<K1*(2*K2+1)+K2;i++){
+  for(i=0;i<USIZE;i++){
    u[i]=u0[i];
  }

@@ -545,69 +593,53 @@ int ns_step_rk4(
  bool keep_en_cst,
  double target_en
 ){
-  int kx,ky;
+  int i;

  // k1
  ns_rhs(tmp1, u, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
  // add to output
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp3[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta/6*tmp1[klookup_sym(kx,ky,K2)];
-    }
+  for(i=0;i<USIZE;i++){
+    tmp3[i]=u[i]+delta/6*tmp1[i];
  }

  // u+h*k1/2
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp2[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta/2*tmp1[klookup_sym(kx,ky,K2)];
-    }
+  for(i=0;i<USIZE;i++){
+    tmp2[i]=u[i]+delta/2*tmp1[i];
  }
  // k2
  ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
  // add to output
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp3[klookup_sym(kx,ky,K2)]+=delta/3*tmp1[klookup_sym(kx,ky,K2)];
-    }
+  for(i=0;i<USIZE;i++){
+    tmp3[i]+=delta/3*tmp1[i];
  }

  // u+h*k2/2
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp2[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta/2*tmp1[klookup_sym(kx,ky,K2)];
-    }
+  for(i=0;i<USIZE;i++){
+    tmp2[i]=u[i]+delta/2*tmp1[i];
  }
  // k3
  ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
  // add to output
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp3[klookup_sym(kx,ky,K2)]+=delta/3*tmp1[klookup_sym(kx,ky,K2)];
-    }
+  for(i=0;i<USIZE;i++){
+    tmp3[i]+=delta/3*tmp1[i];
  }

  // u+h*k3
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp2[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta*tmp1[klookup_sym(kx,ky,K2)];
-    }
+  for(i=0;i<USIZE;i++){
+    tmp2[i]=u[i]+delta*tmp1[i];
  }
  // k4
  ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
  // add to output
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      u[klookup_sym(kx,ky,K2)]=tmp3[klookup_sym(kx,ky,K2)]+delta/6*tmp1[klookup_sym(kx,ky,K2)];
-    }
+  for(i=0;i<USIZE;i++){
+    u[i]=tmp3[i]+delta/6*tmp1[i];
  }

  // keep enstrophy constant
  if(keep_en_cst){
    double en=compute_enstrophy(u, K1, K2, L);
-    for(kx=0;kx<=K1;kx++){
-      for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-	u[klookup_sym(kx,ky,K2)]*=sqrt(target_en/en);
-      }
+    for(i=0;i<USIZE;i++){
+      u[i]*=sqrt(target_en/en);
    }
  }

@@ -634,33 +666,27 @@ int ns_step_rk2(
  bool keep_en_cst,
  double target_en
 ){
-  int kx,ky;
+  int i;

  // k1
  ns_rhs(tmp1, u, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);

  // u+h*k1/2
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp2[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta/2*tmp1[klookup_sym(kx,ky,K2)];
-    }
+  for(i=0;i<USIZE;i++){
+    tmp2[i]=u[i]+delta/2*tmp1[i];
  }
  // k2
  ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
  // add to output
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      u[klookup_sym(kx,ky,K2)]+=delta*tmp1[klookup_sym(kx,ky,K2)];
-    }
+  for(i=0;i<USIZE;i++){
+    u[i]+=delta*tmp1[i];
  }

  // keep enstrophy constant
  if(keep_en_cst){
    double en=compute_enstrophy(u, K1, K2, L);
-    for(kx=0;kx<=K1;kx++){
-      for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-	u[klookup_sym(kx,ky,K2)]*=sqrt(target_en/en);
-      }
+    for(i=0;i<USIZE;i++){
+      u[i]*=sqrt(target_en/en);
    }
  }

@@ -700,7 +726,7 @@ int ns_step_rkf45(
  // whether to compute k1 or leave it as is
  bool compute_k1
 ){
-  int kx,ky;
+  int i;
  double cost;

  // k1: u(t)
@@ -709,53 +735,41 @@ int ns_step_rkf45(
  }

  // k2 : u(t+1/4*delta)
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)/4*k1[klookup_sym(kx,ky,K2)];
-    }
+  for(i=0;i<USIZE;i++){
+    tmp[i]=u[i]+(*delta)/4*k1[i];
  }
  ns_rhs(k2, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);

  // k3 : u(t+3/8*delta)
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(3./32*k1[klookup_sym(kx,ky,K2)]+9./32*k2[klookup_sym(kx,ky,K2)]);
-    }
+  for(i=0;i<USIZE;i++){
+    tmp[i]=u[i]+(*delta)*(3./32*k1[i]+9./32*k2[i]);
  }
  ns_rhs(k3, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);

  // k4 : u(t+12./13*delta)
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(1932./2197*k1[klookup_sym(kx,ky,K2)]-7200./2197*k2[klookup_sym(kx,ky,K2)]+7296./2197*k3[klookup_sym(kx,ky,K2)]);
-    }
+  for(i=0;i<USIZE;i++){
+    tmp[i]=u[i]+(*delta)*(1932./2197*k1[i]-7200./2197*k2[i]+7296./2197*k3[i]);
  }
  ns_rhs(k4, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);

  // k5 : u(t+1*delta)
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(439./216*k1[klookup_sym(kx,ky,K2)]-8*k2[klookup_sym(kx,ky,K2)]+3680./513*k3[klookup_sym(kx,ky,K2)]-845./4104*k4[klookup_sym(kx,ky,K2)]);
-    }
+  for(i=0;i<USIZE;i++){
+    tmp[i]=u[i]+(*delta)*(439./216*k1[i]-8*k2[i]+3680./513*k3[i]-845./4104*k4[i]);
  }
  ns_rhs(k5, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);

  // k6 : u(t+1./2*delta)
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(-8./27*k1[klookup_sym(kx,ky,K2)]+2*k2[klookup_sym(kx,ky,K2)]-3544./2565*k3[klookup_sym(kx,ky,K2)]+1859./4104*k4[klookup_sym(kx,ky,K2)]-11./40*k5[klookup_sym(kx,ky,K2)]);
-    }
+  for(i=0;i<USIZE;i++){
+    tmp[i]=u[i]+(*delta)*(-8./27*k1[i]+2*k2[i]-3544./2565*k3[i]+1859./4104*k4[i]-11./40*k5[i]);
  }
  ns_rhs(k6, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
  
  // next step
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      // u
-      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(25./216*k1[klookup_sym(kx,ky,K2)]+1408./2565*k3[klookup_sym(kx,ky,K2)]+2197./4104*k4[klookup_sym(kx,ky,K2)]-1./5*k5[klookup_sym(kx,ky,K2)]);
-      // U: save to k6, which is no longer needed
-      k6[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(16./135*k1[klookup_sym(kx,ky,K2)]+6656./12825*k3[klookup_sym(kx,ky,K2)]+28561./56430*k4[klookup_sym(kx,ky,K2)]-9./50*k5[klookup_sym(kx,ky,K2)]+2./55*k6[klookup_sym(kx,ky,K2)]);
-    }
+  for(i=0;i<USIZE;i++){
+    // u
+    tmp[i]=u[i]+(*delta)*(25./216*k1[i]+1408./2565*k3[i]+2197./4104*k4[i]-1./5*k5[i]);
+    // U: save to k6, which is no longer needed
+    k6[i]=u[i]+(*delta)*(16./135*k1[i]+6656./12825*k3[i]+28561./56430*k4[i]-9./50*k5[i]+2./55*k6[i]);
  }

  // cost function
@@ -764,10 +778,8 @@ int ns_step_rkf45(
  // compare relative error with tolerance
  if(cost<tolerance){
    // copy to output
-    for(kx=0;kx<=K1;kx++){
-      for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-	u[klookup_sym(kx,ky,K2)]=tmp[klookup_sym(kx,ky,K2)];
-      }
+    for(i=0;i<USIZE;i++){
+      u[i]=tmp[i];
    }
    // next delta to use in future steps (do not exceed max_delta)
    *next_delta=fmin(max_delta, (*delta)*pow(tolerance/cost,0.2));
@@ -775,10 +787,8 @@ int ns_step_rkf45(
    // keep enstrophy constant
    if(keep_en_cst){
      double en=compute_enstrophy(u, K1, K2, L);
-      for(kx=0;kx<=K1;kx++){
-	for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-	  u[klookup_sym(kx,ky,K2)]*=sqrt(target_en/en);
-	}
+      for(i=0;i<USIZE;i++){
+	u[i]*=sqrt(target_en/en);
      }
    }
  }
@@ -825,7 +835,7 @@ int ns_step_rkbs32(
  // whether to compute k1
  bool compute_k1
 ){
-  int kx,ky;
+  int i;
  double cost;

  // k1: u(t)
@@ -835,36 +845,28 @@ int ns_step_rkbs32(
  }

  // k2 : u(t+1/4*delta)
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)/2*(*k1)[klookup_sym(kx,ky,K2)];
-    }
+  for(i=0;i<USIZE;i++){
+    tmp[i]=u[i]+(*delta)/2*(*k1)[i];
  }
  ns_rhs(k2, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);

  // k3 : u(t+3/4*delta)
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(3./4*k2[klookup_sym(kx,ky,K2)]);
-    }
+  for(i=0;i<USIZE;i++){
+    tmp[i]=u[i]+(*delta)*(3./4*k2[i]);
  }
  ns_rhs(k3, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);

  // k4 : u(t+delta)
  // tmp computed here is the next step
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(2./9*(*k1)[klookup_sym(kx,ky,K2)]+1./3*k2[klookup_sym(kx,ky,K2)]+4./9*k3[klookup_sym(kx,ky,K2)]);
-    }
+  for(i=0;i<USIZE;i++){
+    tmp[i]=u[i]+(*delta)*(2./9*(*k1)[i]+1./3*k2[i]+4./9*k3[i]);
  }
  ns_rhs(*k4, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);

  // next step
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      // U: store in k3, which is no longer needed
-      k3[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(7./24*(*k1)[klookup_sym(kx,ky,K2)]+1./4*k2[klookup_sym(kx,ky,K2)]+1./3*k3[klookup_sym(kx,ky,K2)]+1./8*(*k4)[klookup_sym(kx,ky,K2)]);
-    }
+  for(i=0;i<USIZE;i++){
+    // U: store in k3, which is no longer needed
+    k3[i]=u[i]+(*delta)*(7./24*(*k1)[i]+1./4*k2[i]+1./3*k3[i]+1./8*(*k4)[i]);
  }
  
  // compute cost
@@ -873,10 +875,8 @@ int ns_step_rkbs32(
  // compare cost with tolerance
  if(cost<tolerance){
    // add to output
-    for(kx=0;kx<=K1;kx++){
-      for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-	u[klookup_sym(kx,ky,K2)]=tmp[klookup_sym(kx,ky,K2)];
-      }
+    for(i=0;i<USIZE;i++){
+      u[i]=tmp[i];
    }
    // next delta to use in future steps (do not exceed max_delta)
    *next_delta=fmin(max_delta, (*delta)*pow(tolerance/cost,1./3));
@@ -889,10 +889,8 @@ int ns_step_rkbs32(
    // keep enstrophy constant
    if(keep_en_cst){
      double en=compute_enstrophy(u, K1, K2, L);
-      for(kx=0;kx<=K1;kx++){
-	for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-	  u[klookup_sym(kx,ky,K2)]*=sqrt(target_en/en);
-	}
+      for(i=0;i<USIZE;i++){
+	u[i]*=sqrt(target_en/en);
      }
    }
  }
@@ -940,7 +938,7 @@ int ns_step_rkdp54(
  // whether to compute k1
  bool compute_k1
 ){
-  int kx,ky;
+  int i;
  double cost;

  // k1: u(t)
@@ -950,61 +948,47 @@ int ns_step_rkdp54(
  }

  // k2 : u(t+1/5*delta)
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)/5*(*k1)[klookup_sym(kx,ky,K2)];
-    }
+  for(i=0;i<USIZE;i++){
+    tmp[i]=u[i]+(*delta)/5*(*k1)[i];
  }
  ns_rhs(*k2, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);

  // k3 : u(t+3/10*delta)
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(3./40*(*k1)[klookup_sym(kx,ky,K2)]+9./40*(*k2)[klookup_sym(kx,ky,K2)]);
-    }
+  for(i=0;i<USIZE;i++){
+    tmp[i]=u[i]+(*delta)*(3./40*(*k1)[i]+9./40*(*k2)[i]);
  }
  ns_rhs(k3, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);

  // k4 : u(t+4/5*delta)
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(44./45*(*k1)[klookup_sym(kx,ky,K2)]-56./15*(*k2)[klookup_sym(kx,ky,K2)]+32./9*k3[klookup_sym(kx,ky,K2)]);
-    }
+  for(i=0;i<USIZE;i++){
+    tmp[i]=u[i]+(*delta)*(44./45*(*k1)[i]-56./15*(*k2)[i]+32./9*k3[i]);
  }
  ns_rhs(k4, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);

  // k5 : u(t+8/9*delta)
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(19372./6561*(*k1)[klookup_sym(kx,ky,K2)]-25360./2187*(*k2)[klookup_sym(kx,ky,K2)]+64448./6561*k3[klookup_sym(kx,ky,K2)]-212./729*k4[klookup_sym(kx,ky,K2)]);
-    }
+  for(i=0;i<USIZE;i++){
+    tmp[i]=u[i]+(*delta)*(19372./6561*(*k1)[i]-25360./2187*(*k2)[i]+64448./6561*k3[i]-212./729*k4[i]);
  }
  ns_rhs(k5, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);

  // k6 : u(t+delta)
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(9017./3168*(*k1)[klookup_sym(kx,ky,K2)]-355./33*(*k2)[klookup_sym(kx,ky,K2)]+46732./5247*k3[klookup_sym(kx,ky,K2)]+49./176*k4[klookup_sym(kx,ky,K2)]-5103./18656*k5[klookup_sym(kx,ky,K2)]);
-    }
+  for(i=0;i<USIZE;i++){
+    tmp[i]=u[i]+(*delta)*(9017./3168*(*k1)[i]-355./33*(*k2)[i]+46732./5247*k3[i]+49./176*k4[i]-5103./18656*k5[i]);
  }
  ns_rhs(k6, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);

  // k7 : u(t+delta)
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      // tmp computed here is the step
-      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(35./384*(*k1)[klookup_sym(kx,ky,K2)]+500./1113*k3[klookup_sym(kx,ky,K2)]+125./192*k4[klookup_sym(kx,ky,K2)]-2187./6784*k5[klookup_sym(kx,ky,K2)]+11./84*k6[klookup_sym(kx,ky,K2)]);
-    }
+  for(i=0;i<USIZE;i++){
+    // tmp computed here is the step
+    tmp[i]=u[i]+(*delta)*(35./384*(*k1)[i]+500./1113*k3[i]+125./192*k4[i]-2187./6784*k5[i]+11./84*k6[i]);
  }
  // store in k2, which is not needed anymore
  ns_rhs(*k2, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);

  //next step
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      // U: store in k6, which is not needed anymore
-      k6[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(5179./57600*(*k1)[klookup_sym(kx,ky,K2)]+7571./16695*k3[klookup_sym(kx,ky,K2)]+393./640*k4[klookup_sym(kx,ky,K2)]-92097./339200*k5[klookup_sym(kx,ky,K2)]+187./2100*k6[klookup_sym(kx,ky,K2)]+1./40*(*k2)[klookup_sym(kx,ky,K2)]);
-    }
+  for(i=0;i<USIZE;i++){
+    // U: store in k6, which is not needed anymore
+    k6[i]=u[i]+(*delta)*(5179./57600*(*k1)[i]+7571./16695*k3[i]+393./640*k4[i]-92097./339200*k5[i]+187./2100*k6[i]+1./40*(*k2)[i]);
  }
  
  // compute cost
@@ -1013,10 +997,8 @@ int ns_step_rkdp54(
  // compare relative error with tolerance
  if(cost<tolerance){
    // add to output
-    for(kx=0;kx<=K1;kx++){
-      for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-	u[klookup_sym(kx,ky,K2)]=tmp[klookup_sym(kx,ky,K2)];
-      }
+    for(i=0;i<USIZE;i++){
+      u[i]=tmp[i];
    }
    // next delta to use in future steps (do not exceed max_delta)
    *next_delta=fmin(max_delta, (*delta)*pow(tolerance/cost,0.2));
@@ -1029,10 +1011,8 @@ int ns_step_rkdp54(
    // keep enstrophy constant
    if(keep_en_cst){
      double en=compute_enstrophy(u, K1, K2, L);
-      for(kx=0;kx<=K1;kx++){
-	for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-	  u[klookup_sym(kx,ky,K2)]*=sqrt(target_en/en);
-	}
+      for(i=0;i<USIZE;i++){
+	u[i]*=sqrt(target_en/en);
      }
    }
  }
@@ -1137,7 +1117,7 @@ int ns_rhs(
    alpha=compute_alpha(u,K1,K2,g,L);
  }

-  for(i=0; i<K1*(2*K2+1)+K2; i++){
+  for(i=0; i<USIZE; i++){
    out[i]=0;
  }
  for(kx=0;kx<=K1;kx++){
--- a/src/navier-stokes.h
+++ b/src/navier-stokes.h
@@ -1,5 +1,5 @@
 /*
-Copyright 2017-2024 Ian Jauslin
+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.
@@ -31,10 +31,12 @@ typedef struct fft_vects {
 } fft_vect;

 // compute u_k
-int uk( int K1, int K2, int N1, int N2, double final_time, 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, double print_freq, double starting_time, unsigned int nthreadsl, FILE* savefile);
+int uk( int K1, int K2, int N1, int N2, double final_time, 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, double print_freq, double starting_time, unsigned int nthreadsl, FILE* savefile, FILE* utfile, const char* cmd_string, const char* params_string, const char* savefile_string, const char* utfile_string);

 // compute enstrophy and alpha
 int enstrophy( int K1, int K2, int N1, int N2, double final_time, 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, double print_freq, double starting_time, bool print_alpha, unsigned int nthreads, FILE* savefile, FILE* utfile, const char* cmd_string, const char* params_string, const char* savefile_string, const char* utfile_string);
+// compute enstrophy, alpha for the initial condition (useful for debugging)
+int enstrophy_print_init( int K1, int K2, double L, _Complex double* u0, _Complex double* g);

 // compute solution as a function of time, but do not print anything (useful for debugging)
 int quiet( int K1, int K2, int N1, int N2, double final_time, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_cost, double starting_time, _Complex double* u0, _Complex double* g, bool irreversible, bool keep_en_cst, double target_en, unsigned int algorithm, unsigned int nthreads, FILE* savefile);
Author	SHA1	Message	Date
Ian Jauslin	35352a6460	Remove init_enstrophy and init_energy from savefile	2025-12-14 10:29:47 -05:00
Ian Jauslin	8fa9e7f556	Add command to print the enstrophy of the initial condition	2025-10-14 12:20:43 -04:00
Ian Jauslin	afe0498f21	Fix remove_entry	2025-06-05 14:55:34 +02:00
Ian Jauslin	72455cbb65	Print components of u	2025-04-14 18:25:59 -04:00
Ian Jauslin	7471296e59	Typo	2025-03-04 10:57:05 -05:00
Ian Jauslin	08ded444b8	Use savefiles in uk	2025-02-10 15:34:03 -05:00
Ian Jauslin	8a1f3987f4	when reading binary init, only read u unless using the 'resume' command. In particular, removed the 'init_flow' parameter	2025-02-03 14:11:53 -05:00
Ian Jauslin	89791be6d6	Print lyapunov if it is the last step	2025-02-03 12:27:16 -05:00
Ian Jauslin	50c09cf164	Do not sort exponents, and fix norm of flow	2025-02-03 12:08:07 -05:00
Ian Jauslin	21e8dcdb8a	In RK algorithms: do not use kx,ky, use i	2025-02-01 14:12:54 -05:00
Ian Jauslin	e607a4abf9	Ensure that init_flow is used in conjunction with init	2025-02-01 12:15:56 -05:00
Ian Jauslin	6f0f1749a4	Document flow_init in readme	2025-02-01 12:10:50 -05:00
Ian Jauslin	03c2d1b02a	Fix interruption of lyapunov	2025-02-01 11:53:19 -05:00
Ian Jauslin	d1992eca70	interrupting lyapunov mentioned in README	2025-01-31 21:59:41 -05:00
Ian Jauslin	06e5b0e0da	Allow the Lyapunov computation to be interrupted	2025-01-31 21:58:54 -05:00
Ian Jauslin	53a0a0ae4c	Remove D_epsilon	2025-01-31 20:50:55 -05:00
Ian Jauslin	a47ec7896b	Fix bug: every column of tangent flow needs to be evolved	2025-01-31 17:45:25 -05:00
Ian Jauslin	3fa3a86066	Bug fix	2025-01-31 15:04:41 -05:00
Ian Jauslin	0244f03d02	Update copyright	2025-01-31 12:01:53 -05:00
Ian Jauslin	b0f82a2412	fix README layout	2025-01-31 11:09:23 -05:00
Ian Jauslin	170aebfa0c	Lyapunov in README	2025-01-31 11:08:20 -05:00
Ian Jauslin	57df67cc4c	Mention lyapunov algorithms in doc	2025-01-31 10:55:41 -05:00
Ian Jauslin	0f07f025b5	Overhaul of lyapunov exponents	2025-01-31 10:52:40 -05:00
Ian Jauslin	fa52397e87	Correctly deal with complex values in lyapunov exponents documentation	2025-01-27 17:14:46 -05:00
Ian Jauslin	328f3bd99c	New derivation of Lyapunov exponents in doc	2025-01-16 11:42:15 +01:00
Ian Jauslin	acf1c27b1d	Details on RK in doc	2024-12-16 16:57:37 -05:00
Ian Jauslin	9d48b044ac	Reorganize doc	2024-12-15 14:49:30 -05:00