Fix remove_entry

Print components of u
Typo
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 · 2025-02-03 14:11:53 -05:00 · 2025-02-03 12:27:16 -05:00
21 changed files with 2820 additions and 839 deletions
--- a/4
+++ b/4
@ -1,4 +1,4 @@
-# Copyright 2017-2023 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.
@ -30,7 +30,7 @@ LD=$(CC)
 #INCLUDES = 

 #LIBDIRS = 
-LIBS = -lm -lfftw3 -lfftw3_threads
+LIBS = -lm -lfftw3 -lfftw3_threads -lopenblas_64


 override LDFLAGS +=$(LIBDIRS)$(LIBS)
--- a/README.md
+++ b/README.md
@ -39,9 +39,22 @@ The available commands are

 * `enstrophy`: to compute the enstrophy and various other observables. This
  command prints
-  ```step_index time average(alpha) average(enstrophy) average(alpha*enstrophy) alpha enstrophy alpha*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. 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.

@ -100,8 +113,13 @@ should be a `;` sperated list of `key=value` pairs. The possible keys are
  force, excluding kx<0 and kx=0&&ky<=0. The plaintext file format is
  `kx ky real_part imag_part`.

-* `init_en` (double, default 1.54511597324389e+02): initial value of the energy if
-  `equation=irreversible` or of the enstrophy if `equation=reversible`.
+* `init_energy` (double, default is to not rescale the initial condition, is
+  incompatible with `init_enstrophy`): enforce a value for the energy of the
+  initial condition.
+
+* `init_enstrophy` (double, default is to not rescale the initial condition, is
+  incompatible with `init_energy`): enforce a value for the enstrophy of the
+  initial condition.

 * `random_seed` (int, default ): seed for random initialization.

@ -120,28 +138,46 @@ should be a `;` sperated list of `key=value` pairs. The possible keys are
 * `max_delta` (double, default 1e-2): when using the adaptive step, do not
  exceet `delta_max`.

-* `adaptive_norm`: norm to use to estimate the error of the adaptive method:
-  `L1` (default) for the normalized L1 norm, `k3` for the normalized L1 norm of
-  f_k/|k|^3, `k32` for the normalized L1 norm, `enstrophy` for the enstrophy
-  norm.
+* `adaptive_cost`: cost function to use to estimate the error of the adaptive
+  method: `L1` (default) for the normalized L1 norm, `k3` for the normalized L1
+  norm of f_k/|k|^3, `k32` for the normalized L1 norm, `enstrophy` for the
+  enstrophy, `alpha` for alpha.

 * `keep_en_cst` (0 or 1, default 0): impose that the enstrophy is constant at
  each step (only really useful for the reversible equation).

+* `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-2023 Ian Jauslin
+Copyright 2017-2025 Ian Jauslin
--- a/docs/nstrophy_doc.tex
+++ b/docs/nstrophy_doc.tex
@ -1,4 +1,4 @@
-% Copyright 2017-2023 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.
@ -15,6 +15,7 @@
 \documentclass{ian}

 \usepackage{largearray}
+\usepackage{dsfont}

 \begin{document}

@ -86,6 +87,7 @@ and $q\cdot p^\perp$ is antisymmetric under exchange of $q$ and $p$. Therefore,
  \partial_t\hat u_k=
  -\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)
+  =:\mathfrak F_k(\hat u)
  \label{ins_k}
 \end{equation}
 with
@ -103,104 +105,7 @@ We truncate the Fourier modes and assume that $\hat u_k=0$ if $|k_1|>K_1$ or $|k
 \end{equation}
 \bigskip

-\point{\bf 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)$.
-{\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}
-  \item the Runge-Kutta-Dormand-Prince method ({\tt RKDP54}), which is of 5th order, and adjusts the step by comparing to a 4th order method;
-  \item the Runge-Kutta-Fehlberg method ({\tt RKF45}), which is of 4th order, and adjusts the step by comparing to a 5th order method;
-  \item the Runge-Kutta-Bogacki-Shampine method ({\tt RKBS32}), which is of 3d order, and adjusts the step by comparing to a 2nd order method.
-\end{itemize}
-In these adaptive step methods, two steps are computed at different orders: $\hat u_k^{(n)}$ and $\hat U_k^{(n)}$, the step size is adjusted at every step in such a way that the error is small enough:
-\begin{equation}
-  \|\hat u^{(n)}-\hat U^{(n)}\|
-  <\epsilon_{\mathrm{target}}
-\end{equation}
-for some given $\epsilon_{\mathrm{target}}$, set using the {\tt adaptive\_tolerance} parameter.
-The choice of the norm matters, and will be discussed below.
-If the error is larger than the target, then the step size is decreased.
-How this is done depends on the order of algorithm.
-If the order is $q$ (here we mean the smaller of the two orders, so 4 for {\tt RKDP54} and {\tt RKF45} and 2 for {\tt RKBS32}), then we expect
-\begin{equation}
-  \|\hat u^{(n)}-\hat U^{(n)}\|=\delta_n^qC_n
-  .
-\end{equation}
-We wish to set $\delta_{n+1}$ so that
-\begin{equation}
-  \delta_{n+1}^qC_n=\epsilon_{\mathrm{target}}
-\end{equation}
-so
-\begin{equation}
-  \delta_{n+1}
-  =\left(\frac{\epsilon_{\mathrm{target}}}{C_n}\right)^{\frac1q}
-  =\delta_n\left(\frac{\epsilon_{\mathrm{target}}}{\|\hat u^{(n)}-\hat U^{(n)}\|}\right)^{\frac1q}
-  .
-  \label{adaptive_delta}
-\end{equation}
-(Actually, to be safe and ensure that $\delta$ decreases sufficiently, we multiply this by a safety factor that can be set using the {\tt adaptive\_factor} parameter.)
-If the error is smaller than the target, we increase $\delta$ using\-~(\ref{adaptive_delta}) (without the safety factor).
-To be safe, we also set a maximal value for $\delta$ via the {\tt max\_delta} parameter.
-\bigskip
-
-\indent
-The choice of the norm $\|\cdot\|$ matters.
-It can be made by specifying the parameter {\tt adaptive\_norm}.
-\begin{itemize}
-  \item A naive choice is to take $\|\cdot\|$ to be the normalized $L_1$ norm:
-  \begin{equation}
-    \|f\|:=
-    \frac1{\mathcal N}\sum_k|f_k|
-    ,\quad
-    \mathcal N:=\sum_k|\hat u_k^{(n)}-\hat u_k^{(n-1)}|
-    .
-  \end{equation}
-  This norm is selected by choosing {\tt adaptive\_norm=L1}.
-
-  \item Empirically, we have found that $|\hat u-\hat U|$ behaves like $k^{-3}$ for {\tt RKDP54} and {\tt RKF45}, and like $k^{-\frac32}$ for {\tt RKBS32}, so a norm of the form
-  \begin{equation}
-    \|f\|:=\frac1{\mathcal N}\sum_k|f_k|k^{-3}
-    ,\quad
-    \mathcal N:=\sum_k|\hat u_k^{(n)}-\hat u_k^{(n-1)}|k^{-3}
-  \end{equation}
-  or
-  \begin{equation}
-    \|f\|:=\frac1{\mathcal N}\sum_k|f_k|k^{-\frac32}
-    ,\quad
-    \mathcal N:=\sum_k|\hat u_k^{(n)}-\hat u_k^{(n-1)}|k^{-\frac32}
-  \end{equation}
-  are sensible choices.
-  These norms are selected by choosing {\tt adaptive\_norm=k3} and {\tt adaptive\_norm=k32} respectively.
-
-  \item
-  Another option is to define a norm based on the expression of the enstrophy\-~(\ref{enstrophy}):
-  \begin{equation}
-    \|f\|:=\frac1{\mathcal N}\sqrt{\sum_k k^2|f_k|^2}
-    ,\quad
-    \mathcal N:=\frac{\sqrt{\sum_k k^2|\hat u_k^{(n)}|^2}+\sqrt{\sum_k k^2|\hat U_k^{(n)}|^2}}{\sum_k k^2|\hat u_k^{(n)}|^2}
-    .
-  \end{equation}
-  Doing so controls the error of the enstrophy through
-  \begin{equation}
-    \frac1{\mathcal N^2}|\mathcal En(\hat u)-\mathcal En(\hat U)|\equiv|\|\hat u\|^2-\|\hat U\|^2|\leqslant \|\hat u-\hat U\|(\|\hat u\|+\|\hat U\|)
-  \end{equation}
-  so
-  \begin{equation}
-    \frac1{\mathcal N^2}
-    |\mathcal En(\hat u)-\mathcal En(\hat U)|\leqslant
-    \|\hat u-\hat U\|\frac1{\mathcal N}\left(\sqrt{\sum_k k^2|\hat u_k|^2}+\sqrt{\sum_k k^2|\hat U_k|^2}\right)
-  \end{equation}
-  and thus
-  \begin{equation}
-    \frac{|\mathcal En(\hat u)-\mathcal En(\hat U)|}{\mathcal En(\hat u)}\leqslant
-    \|\hat u-\hat U\|
-    .
-  \end{equation}
-  This norm is selected by choosing {\tt adaptive\_norm=enstrophy}.
-\end{itemize}
-\bigskip
-
-\point{\bf Reality}.
+{\bf Remark}:
 Since $U$ is real, $\hat U_{-k}=\hat U_k^*$, and so
 \begin{equation}
  \hat u_{-k}=\hat u_k^*
@ -222,9 +127,159 @@ Thus,
  \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

-\point{\bf FFT}. We compute T using a fast Fourier transform, defined as
+\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}
+  \item the Runge-Kutta-Dormand-Prince method ({\tt RKDP54}), which is of 5th order, and adjusts the step by comparing to a 4th order method;
+  \item the Runge-Kutta-Fehlberg method ({\tt RKF45}), which is of 4th order, and adjusts the step by comparing to a 5th order method;
+  \item the Runge-Kutta-Bogacki-Shampine method ({\tt RKBS32}), which is of 3d order, and adjusts the step by comparing to a 2nd order method.
+\end{itemize}
+In these adaptive step methods, two steps are computed at different orders: $\hat u_k^{(n)}$ and $\hat U_k^{(n)}$, the step size is adjusted at every step in such a way that the error is small enough:
+\begin{equation}
+  D(\hat u^{(n)},\hat U^{(n)})
+  <\epsilon_{\mathrm{target}}
+\end{equation}
+for some given {\it cost function} $D$, and $\epsilon_{\mathrm{target}}$, set using the {\tt adaptive\_tolerance} parameter.
+The choice of the cost function matters, and will be discussed below.
+If the error is larger than the target, then the step size is decreased.
+How this is done depends on the order of algorithm.
+If the order is $q$ (here we mean the smaller of the two orders, so 4 for {\tt RKDP54} and {\tt RKF45} and 2 for {\tt RKBS32}), then we expect (as long as the cost function is such that $D(\hat u,\hat u+\varphi)\sim\|\varphi\|$ in some norm)
+\begin{equation}
+  D(\hat u^{(n)},\hat U^{(n)})=\delta_n^qC_n
+\end{equation}
+for some number $C_n$.
+We wish to set $\delta_{n+1}$ so that
+\begin{equation}
+  \delta_{n+1}^qC_n=\epsilon_{\mathrm{target}}
+\end{equation}
+so
+\begin{equation}
+  \delta_{n+1}
+  =\left(\frac{\epsilon_{\mathrm{target}}}{C_n}\right)^{\frac1q}
+  =\delta_n\left(\frac{\epsilon_{\mathrm{target}}}{D(\hat u^{(n)},\hat U^{(n)})}\right)^{\frac1q}
+  .
+  \label{adaptive_delta}
+\end{equation}
+(Actually, to be safe and ensure that $\delta$ decreases sufficiently, we multiply this by a safety factor that can be set using the {\tt adaptive\_factor} parameter.)
+If the error is smaller than the target, we increase $\delta$ using\-~(\ref{adaptive_delta}) (without the safety factor).
+To be safe, we also set a maximal value for $\delta$ via the {\tt max\_delta} parameter.
+\bigskip
+
+\indent
+The choice of the cost function $D$ matters.
+It can be made by specifying the parameter {\tt adaptive\_cost}.
+\begin{itemize}
+  \item
+  For computations where the main focus is the enstrophy\-~(\ref{enstrophy}), one may want to set the cost function to the relative difference of the enstrophies:
+  \begin{equation}
+    D(\hat u,\hat U):=\frac{|\mathcal En(\hat u)-\mathcal En(\hat U)|}{\mathcal En(\hat u)}
+    .
+  \end{equation}
+  This cost function is selected by choosing {\tt adaptive\_cost=enstrophy}.
+
+  \item
+  For computations where the main focus is the value of $\alpha$\-~(\ref{alpha}), one may want to set the cost function to the relative difference of $\alpha$:
+  \begin{equation}
+    D(\hat u,\hat U):=\frac{|\alpha(\hat u)-\alpha(\hat U)|}{|\alpha(\hat u)|}
+    .
+  \end{equation}
+  This cost function is selected by choosing {\tt adaptive\_cost=alpha}.
+
+  \item Alternatively, one my take $D$ to be the normalized $L_1$ norm:
+  \begin{equation}
+    D(\hat u,\hat U):=
+    \frac1{\mathcal N}\sum_k|\hat u_k-\hat U_k|
+    ,\quad
+    \mathcal N:=\sum_k|\hat u_k|
+    .
+  \end{equation}
+  This function is selected by choosing {\tt adaptive\_cost=L1}.
+
+  \item Empirically, we have found that $|\hat u-\hat U|$ behaves like $k^{-3}$ for {\tt RKDP54} and {\tt RKF45}, and like $k^{-\frac32}$ for {\tt RKBS32}, so a cost function of the form
+  \begin{equation}
+    D(\hat u,\hat U):=\frac1{\mathcal N}\sum_k|\hat u_k-\hat U_k|k^{-3}
+    ,\quad
+    \mathcal N:=\sum_k|\hat u_k|k^{-3}
+  \end{equation}
+  or
+  \begin{equation}
+    D(\hat u,\hat U):=\frac1{\mathcal N}\sum_k|\hat u_k-\hat U_k|k^{-\frac32}
+    ,\quad
+    \mathcal N:=\sum_k|\hat u_k|k^{-\frac32}
+  \end{equation}
+  are sensible choices.
+  These cost functions are selected by choosing {\tt adaptive\_cost=k3} and {\tt adaptive\_cost=k32} respectively.
+\end{itemize}
+
+
+\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}
@ -246,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}
@ -267,9 +322,13 @@ Therefore,
    \mathcal F\left(q_x|q|\hat u_q\right)(n)
  \right)(k)
 \end{equation}
+
+\subsection{Observables}
+\indent
+We define the following observables.
 \bigskip

-\point{\bf Energy}.
+\subsubsection{Energy}.
 We define the energy as
 \begin{equation}
  E(t)=\frac12\int\frac{dx}{L^2}\ U^2(t,x)=\frac12\sum_{k\in\mathbb Z^2}|\hat U_k|^2
@ -370,73 +429,243 @@ and
  .
  \label{enstrophy}
 \end{equation}
-\bigskip

-\point{\bf Enstrophy}.
+\subsubsection{Enstrophy}.
 The enstrophy is defined as
 \begin{equation}
  \mathcal En(t)=\int\frac{dx}{L^2}\ |\nabla U|^2
  =\frac{4\pi^2}{L^2}\sum_{k\in\mathbb Z^2}k^2|\hat U_k|^2
  .
 \end{equation}
+
+\subsection{Lyapunov exponents}
+\indent
+The Lyapunov are defined from the {\it tangent flow} of the dynamics.
+Consider an equation of the form
+\begin{equation}
+  \dot u=f(t;u)
+  .
+\end{equation}
+Now, the flow may not be complex-differentiable, so the tangent flow should be computed on the real and imaginary parts.
+Let
+\begin{equation}
+  u=\zeta+i\xi
+  ,\quad
+  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

-\point{\bf 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:
+\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}
-  \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
+  [0,t)=\bigcup_{i=0}^{N-1}[t_i,t_{i+1})
  ,\quad
-  k\cdot\hat U_k=0
+  \varphi_{0,t}=\prod_{i=N-1}^0\varphi_{t_{i},t_{i+1}}
+  .
 \end{equation}
-where $\alpha$ is chosen such that the enstrophy is constant.
-In terms of $\hat u$\-~(\ref{udef}), (\ref{gdef}), (\ref{T}):
+At the thresholds between these intervals, we perform a QR decomposition:
 \begin{equation}
-  \partial_t\hat u_k=
+  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}
+%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}
+For the reversible equation,
+\begin{equation}
+  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}
-  .
-\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-2023 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-2023 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-2023 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.
@ -14,10 +14,13 @@ See the License for the specific language governing permissions and
 limitations under the License.
 */

+#define M_PI 3.14159265358979323846

 #define COMMAND_UK 1
 #define COMMAND_ENSTROPHY 2
 #define COMMAND_QUIET 3
+#define COMMAND_RESUME 4
+#define COMMAND_LYAPUNOV 5

 #define DRIVING_ZERO 1
 #define DRIVING_TEST 2
@ -37,7 +40,14 @@ limitations under the License.
 #define ALGORITHM_RKDP54 1002
 #define ALGORITHM_RKBS32 1003

-#define NORM_L1 1
-#define NORM_k3 2
-#define NORM_k32 3
-#define NORM_ENSTROPHY 4
+#define COST_L1 1
+#define COST_k3 2
+#define COST_k32 3
+#define COST_ENSTROPHY 4
+#define COST_ALPHA 5
+
+#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-2023 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-2023 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
@ -0,0 +1,244 @@
+/*
+Copyright 2017-2025 Ian Jauslin
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+#include "dstring.h"
+#include <stdarg.h>
+#include <stdlib.h>
+#include <string.h>
+
+// init
+int dstring_init (dstring* str, unsigned int memory){
+  str->string=calloc(memory,sizeof(char));
+  str->memory=memory;
+  str->length=0;
+  return(0);
+}
+int dstring_free (dstring str){
+  free(str.string);
+  return(0);
+}
+
+// resize memory
+int dstring_resize (dstring* str, unsigned int newsize){
+  unsigned int i;
+  dstring new_str;
+  dstring_init (&new_str, newsize);
+  for(i=0;i<str->length;i++){
+    new_str.string[i]=str->string[i];
+  }
+  new_str.length=str->length;
+  free(str->string);
+  *str=new_str;
+  return(0);
+}
+
+// add a character
+int dstring_append (char val, dstring* output){
+  if(output->length >= output->memory){
+    dstring_resize (output,2*output->memory+1);
+  }
+  output->string[output->length]=val;
+  output->length++;
+  return(0);
+}
+
+// copy
+int dstring_cpy (dstring input, dstring* output){
+  dstring_init (output, input.length);
+  dstring_cpy_noinit (input, output);
+  return(0);
+}
+// do not init output str
+int dstring_cpy_noinit (dstring input, dstring* output){
+  unsigned int i;
+  if(output->memory<input.length){
+    fprintf(stderr,"error: cannot copy a dstring of length %u to a dstring with memory %u\n",input.length,output->memory);
+    return(-1);
+  }
+  for(i=0;i<input.length;i++){
+    output->string[i]=input.string[i];
+  }
+  output->length=input.length;
+  return(0);
+}
+
+
+// concatenate
+int dstring_concat (dstring input, dstring* output){
+  unsigned int i;
+  for(i=0;i<input.length;i++){
+    dstring_append (input.string[i], output);
+  }
+  return(0);
+}
+
+// sub-str
+int dstring_substr (dstring str, unsigned int begin, unsigned int end, dstring* substr){
+  unsigned int i;
+  if(begin>end || end>=str.length){
+    fprintf(stderr,"error: cannot take substring [%u,%u] of dstring of lengdth %u\n",begin,end,str.length);
+    return(-1);
+  }
+  dstring_init (substr,end-begin);
+  for(i=begin;i<=end;i++){
+    dstring_append (str.string[i], substr);
+  }
+  return(0);
+}
+
+// find
+int dstring_find (char val, dstring str){
+  unsigned int i;
+  for(i=0;i<str.length;i++){
+    if(str.string[i]==val){
+      return(i);
+    }
+  }
+  return(-1);
+}
+
+// compare strings
+int dstring_cmp (dstring str1, dstring str2){
+  unsigned int i;
+
+  // compare lengths
+  if(str1.length<str2.length){
+    return(-1);
+  }
+  if(str1.length>str2.length){
+    return(1);
+  }
+
+  // compare terms
+  for(i=0;i<str1.length;i++){
+    if(str1.string[i]<str2.string[i]){
+      return(-1);
+    }
+    if(str1.string[i]<str2.string[i]){
+      return(1);
+    }
+  }
+
+  // if equal
+  return(0);
+}
+
+// print str
+int dstring_print (dstring str, FILE* file){
+  unsigned int i;
+  for(i=0;i<str.length;i++){
+    fprintf(file, "%s",str.string);
+  }
+  return(0);
+}
+
+// append a char*
+int dstring_append_str(char* str, dstring* output){
+  char* ptr;
+  for (ptr=str;*ptr!='\0';ptr++){
+    dstring_append(*ptr, output);
+  }
+  return(0);
+}
+
+// convert to char*
+int dstring_to_str(dstring input, char** output){
+  unsigned int i;
+  (*output)=calloc(input.length+1,sizeof(char));
+  for(i=0;i<input.length;i++){
+    (*output)[i]=input.string[i];
+  }
+  if((*output)[input.length-1]!='\0'){
+    (*output)[input.length]='\0';
+  }
+  return(0);
+}
+// noinit (changes the size of input if needed)
+char* dstring_to_str_noinit(dstring* input){
+  // check if string ends in a trailing '\0' (or if string is empty)
+  if(input->length==0 || input->string[input->length-1]!='\0'){
+    if(input->length==input->memory){
+      dstring_resize(input,input->length+1);
+    }
+    // add final '\0'
+    input->string[input->length]='\0';
+  }
+  return(input->string);
+}
+
+// convert from char*
+int str_to_dstring(char* str, dstring* output){
+  dstring_init(output, strlen(str));
+  str_to_dstring_noinit(str, output);
+  return(0);
+}
+int str_to_dstring_noinit(char* str, dstring* output){
+  // reset length
+  output->length=0;
+  dstring_append_str(str, output);
+  return(0);
+}
+
+// compare a dstring and a char*
+int dstring_cmp_str(dstring dstring, char* str){
+  unsigned int j;
+  for(j=0;j<dstring.length && str[j]!='\0';j++){
+    if(dstring.string[j]!=str[j]){
+      return(1);
+    }
+  }
+  if(j==dstring.length && str[j]=='\0'){
+    return(0);
+  }
+  return(1);
+}
+
+// format strings
+int dstring_snprintf(dstring* output, char* fmt, ...){
+  size_t size=100;
+  unsigned int extra_size;
+  char* out_str=calloc(size,sizeof(char));
+  char* ptr;
+  va_list vaptr;
+
+  // initialize argument list starting after fmt
+  va_start(vaptr, fmt);
+  // print format
+  extra_size=vsnprintf(out_str, size, fmt, vaptr);
+  va_end(vaptr);
+
+  // if too large
+  if(extra_size>size){
+    // resize
+    free(out_str);
+    // +1 for '\0'
+    size=extra_size+1;
+    out_str=calloc(size,sizeof(char));
+    // read format again
+    va_start(vaptr, fmt);
+    vsnprintf(out_str,size,fmt,vaptr);
+    va_end(vaptr);
+  }
+
+  // write to char array
+  for(ptr=out_str;*ptr!='\0';ptr++){
+    dstring_append(*ptr, output);
+  }
+
+  free(out_str);
+
+  return(0);
+}
--- a/src/dstring.h
+++ b/src/dstring.h
@ -0,0 +1,77 @@
+/*
+Copyright 2017-2025 Ian Jauslin
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+#ifndef DSTRING_H
+#define DSTRING_H
+
+#include <stdio.h>
+
+typedef struct dstring {
+  char* string;
+  unsigned int length;
+  unsigned int memory;
+} dstring;
+
+// init
+int dstring_init (dstring* str, unsigned int memory);
+int dstring_free (dstring str);
+
+// resize memory
+int dstring_resize (dstring* str, unsigned int newsize);
+
+// add a character
+int dstring_append (char val, dstring* output);
+
+// copy
+int dstring_cpy (dstring input, dstring* output);
+// do not init output str
+int dstring_cpy_noinit (dstring input, dstring* output);
+
+
+// concatenate
+int dstring_concat (dstring input, dstring* output);
+
+// sub-str
+int dstring_substr (dstring str, unsigned int begin, unsigned int end, dstring* substr);
+
+// find
+int dstring_find (char val, dstring str);
+
+// compare strings
+int dstring_cmp (dstring str1, dstring str2);
+
+// print str
+int dstring_print (dstring str, FILE* file);
+
+// append a char*
+int dstring_append_str(char* str, dstring* output);
+
+// convert to char*
+int dstring_to_str(dstring input, char** output);
+// noinit (changes the size of input if needed)
+char* dstring_to_str_noinit(dstring* input);
+
+// convert from char*
+int str_to_dstring(char* str, dstring* output);
+int str_to_dstring_noinit(char* str, dstring* output);
+
+// compare a dstring and a char*
+int dstring_cmp_str(dstring dstring, char* str);
+
+// format strings
+int dstring_snprintf(dstring* output, char* fmt, ...);
+
+#endif
--- a/src/init.c
+++ b/src/init.c
@ -1,5 +1,5 @@
 /*
-Copyright 2017-2023 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.
@ -17,21 +17,17 @@ limitations under the License.
 #include "init.h"
 #include "navier-stokes.h"
 #include "io.h"
-#include <stdlib.h>
 #include <math.h>
+#include <stdlib.h>

 // random initial condition
 int init_random (
  _Complex double* u0,
-  double init_en,
  int K1,
  int K2,
-  double L,
-  int seed,
-  bool irreversible
+  int seed
 ){
  int kx,ky;
-  double rescale;
  double x,y;

  srand(seed);
@ -45,33 +41,16 @@ int init_random (
    }
  }

-  if (irreversible) {
-    // fix the energy
-    rescale=compute_energy(u0, K1, K2);
-  } else {
-    // fix the enstrophy
-    rescale=compute_enstrophy(u0, K1, K2, L);
-  }
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      u0[klookup_sym(kx,ky,K2)]=u0[klookup_sym(kx,ky,K2)]*sqrt(init_en/rescale);
-    }
-  }
-
  return 0;
 }

 // Gaussian initial condition
 int init_gaussian (
  _Complex double* u0,
-  double init_en,
  int K1,
-  int K2,
-  double L,
-  bool irreversible
+  int K2
 ){
  int kx,ky;
-  double rescale;

  for(kx=0;kx<=K1;kx++){
    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
@ -79,19 +58,6 @@ int init_gaussian (
    }
  }

-  if (irreversible) {
-    // fix the energy
-    rescale=compute_energy(u0, K1, K2);
-  } else {
-    // fix the enstrophy
-    rescale=compute_enstrophy(u0, K1, K2, L);
-  }
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      u0[klookup_sym(kx,ky,K2)]=u0[klookup_sym(kx,ky,K2)]*sqrt(init_en/rescale);
-    }
-  }
-
  return 0;
 }

--- a/src/init.h
+++ b/src/init.h
@ -1,5 +1,5 @@
 /*
-Copyright 2017-2023 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,10 +21,10 @@ limitations under the License.
 #include <stdbool.h>

 // random initial condition
-int init_random(_Complex double* u0, double init_en, int K1, int K2, double L, int seed, bool irreversible);
+int init_random(_Complex double* u0, int K1, int K2, int seed);

 // Gaussian initial condition
-int init_gaussian(_Complex double* u0, double init_en, int K1, int K2, double L, bool irreversible);
+int init_gaussian(_Complex double* u0, int K1, int K2);

 // Initialize from file
 int init_file_txt (_Complex double* u0, int K1, int K2, FILE* initfile);
--- a/src/int_tools.c
+++ b/src/int_tools.c
@ -1,5 +1,5 @@
 /*
-Copyright 2017-2023 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-2023 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-2023 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.
@ -17,6 +17,7 @@ limitations under the License.
 #include <stdbool.h>
 #include <stdlib.h>
 #include <string.h>
+#include "constants.cpp"
 #include "io.h"
 #include "navier-stokes.h"

@ -50,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;
@ -70,7 +95,7 @@ int read_vec(_Complex double* out, int K1, int K2, FILE* file){
  }

  // allocate line buffer
-  line=calloc(sizeof(char), len);
+  line=calloc(len,sizeof(char));

  while(1){
    c=fgetc(file);
@ -125,7 +150,7 @@ int read_vec(_Complex double* out, int K1, int K2, FILE* file){
      // check that there is room in buffer
      if(pos==len){
 	// too short: reallocate
-	line_realloc=calloc(sizeof(char), 2*len);
+	line_realloc=calloc(2*len,sizeof(char));
 	for(pos=0;pos<len;pos++){
 	  line_realloc[pos]=line[pos];
 	}
@ -148,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=='#'){
@ -183,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;
 }

@ -197,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';
@ -225,3 +301,119 @@ int remove_entry(

  return 0;
 }
+
+
+// save state to savefile
+int save_state(
+  _Complex double* u,
+  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
+){
+  if(savefile==NULL){
+    fprintf(stderr, "error while writing savefile: savefile not open\n");
+    return -1;
+  }
+
+  fprintf(savefile,"# Continue computation with\n");
+
+  // command to resume
+  fprintf(savefile,"#! ");
+  fprintf(savefile, cmd_string);
+  // params
+  // allocate buffer for params
+  if(params_string!=NULL) {
+    char* params=calloc(strlen(params_string)+1,sizeof(char));
+    strcpy(params, params_string);
+    remove_entry(params, "starting_time");
+    remove_entry(params, "init");
+    if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){
+      remove_entry(params, "delta");
+    }
+    fprintf(savefile," -p \"%s;init=file:%s", params, savefile_string);
+    free(params);
+    // write delta if adaptive, and not writing binary
+    if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD && (savefile==stderr || savefile==stdout)){
+      fprintf(savefile,";delta=%.15e", step);
+    }
+    // write starting_time if not writing binary
+    if(savefile==stderr || savefile==stdout){
+      fprintf(savefile,";starting_time=%.15e", time);
+    }
+    fprintf(savefile,"\"");
+  }
+
+  // utfile
+  if(utfile!=NULL){
+    fprintf(savefile," -u \"%s\"", utfile_string);
+  }
+  // threads
+  if(nthreads!=1){
+    fprintf(savefile," -t %d", nthreads);
+  }
+
+  switch (command){
+  case COMMAND_UK:
+    fprintf(savefile," uk\n");
+    break;
+  case COMMAND_ENSTROPHY:
+    fprintf(savefile," enstrophy\n");
+    break;
+  case COMMAND_QUIET:
+    fprintf(savefile," quiet\n");
+    break;
+  case COMMAND_LYAPUNOV:
+    fprintf(savefile," lyapunov\n");
+    break;
+  case COMMAND_RESUME:
+    fprintf(savefile," resume\n");
+    break;
+
+  default:
+    fprintf(stderr,"error while writing savefile: unrecognized command %d\n", command);
+    return(-1);
+    break;
+  }
+
+  // save final u to savefile
+  if(savefile==stderr || savefile==stdout){
+    write_vec(u, K1, K2, savefile);
+  } else {
+    write_vec_bin(u, K1, K2, savefile);
+    // last binary entry: starting time
+    fwrite(&time, sizeof(double), 1, savefile);
+    // extra binary data for adaptive algorithm
+    if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){
+      fwrite(&step, sizeof(double), 1, savefile);
+    }
+  }
+
+  return 0;
+}
+
+// write u to utfile in plain text
+int write_utfile(
+  _Complex double* u,
+  FILE* utfile,
+  int K1,
+  int K2
+){
+
+  if(utfile==NULL){
+    fprintf(stderr, "error while writing utfile: utfile is not open\n");
+    return -1;
+  }
+
+  write_vec(u, K1, K2, utfile);
+  return 0;
+}
--- a/src/io.h
+++ b/src/io.h
@ -1,5 +1,5 @@
 /*
-Copyright 2017-2023 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,12 +22,28 @@ 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);

+// save state to savefile
+int save_state(_Complex double* u, 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);
+
+// write u to utfile in plain text
+int write_utfile(_Complex double* u, FILE* utfile, int K1, int K2);
 #endif
--- a/src/lyapunov.c
+++ b/src/lyapunov.c
@ -0,0 +1,815 @@
+/*
+Copyright 2017-2025 Ian Jauslin
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+#include "constants.cpp"
+#include "lyapunov.h"
+#include "io.h"
+#include <openblas64/cblas.h>
+#include <openblas64/lapacke.h>
+#include <math.h>
+#include <stdlib.h>
+
+#define MATSIZE (2*(K1*(2*K2+1)+K2))
+#define USIZE (K1*(2*K2+1)+K2)
+
+// compute Lyapunov exponents
+int lyapunov(
+  int K1,
+  int K2,
+  int N1,
+  int N2,
+  double final_time,
+  double lyapunov_reset,
+  unsigned int lyapunov_trigger,
+  double nu,
+  double delta,
+  double L,
+  double adaptive_tolerance,
+  double adaptive_factor,
+  double max_delta,
+  unsigned int adaptive_cost,
+  _Complex double* u0,
+  _Complex double* g,
+  bool irreversible,
+  bool keep_en_cst,
+  double target_en,
+  unsigned int algorithm,
+  unsigned int algorithm_lyapunov,
+  double starting_time,
+  unsigned int nthreads,
+  double* flow0,
+  double* lyapunov_avg0,
+  double prevtime, // the previous time at which a QR decomposition was performed
+  double lyapunov_startingtime, // the time at which the lyapunov exponent computation was started (useful in interrupted computation)
+  FILE* savefile,
+  FILE* utfile,
+  // for interrupt recovery
+  const char* cmd_string,
+  const char* params_string,
+  const char* savefile_string,
+  const char* utfile_string
+){
+  double* lyapunov;
+  double* lyapunov_avg;
+  double* flow;
+  _Complex double* u;
+  double time;
+  fft_vect fftu1;
+  fft_vect fftu2;
+  fft_vect fftu3;
+  fft_vect fftu4;
+  fft_vect fft1;
+  fft_vect ifft;
+  double* tmp1;
+  double* tmp2;
+  double* tmp3;
+  _Complex double* tmp4;
+  _Complex double* tmp5;
+  _Complex double* tmp6;
+  _Complex double* tmp7;
+  _Complex double* tmp8;
+  _Complex double* tmp9;
+  _Complex double* tmp10;
+  double* tmp11;
+  int i,j;
+
+  // period (only useful with LYAPUNOV_TRIGGER_TIME, but compute it anyways: it won't take much time...)
+  // add 0.1 to ensure proper rounding
+  uint64_t n=(uint64_t)((starting_time-fmod(starting_time, lyapunov_reset))/lyapunov_reset+0.1);
+
+  lyapunov_init_tmps(&lyapunov, &lyapunov_avg, &flow, &u, &tmp1, &tmp2, &tmp3, &tmp4, &tmp5, &tmp6, &tmp7, &tmp8, &tmp9, &tmp10, &tmp11, &fftu1, &fftu2, &fftu3, &fftu4, &fft1, &ifft, K1, K2, N1, N2, nthreads, algorithm, algorithm_lyapunov);
+
+  // copy initial condition
+  copy_u(u, u0, K1, K2);
+
+  // initialize flow and averages
+  for (i=0;i<MATSIZE;i++){
+    for (j=0;j<MATSIZE;j++){
+      flow[i*MATSIZE+j]=flow0[i*MATSIZE+j];
+    }
+    lyapunov_avg[i]=lyapunov_avg0[i];
+  }
+
+  // store step (useful for adaptive step methods)
+  double step=delta;
+  double next_step=step;
+
+  // iterate
+  time=starting_time;
+  while(final_time<0 || time<final_time){
+    // 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;
+
+    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_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(flow[i*MATSIZE+i]))/(time-prevtime);
+      }
+
+      //// 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>lyapunov_startingtime)){
+	  lyapunov_avg[i]=lyapunov_avg[i]*(prevtime-lyapunov_startingtime)/(time-lyapunov_startingtime)+lyapunov[i]*(time-prevtime)/(time-lyapunov_startingtime);
+	}
+      }
+
+      // print
+      for(i=0; i<MATSIZE; i++){
+	printf("% .15e % .15e % .15e\n",time, lyapunov[i], lyapunov_avg[i]);
+      }
+      printf("\n");
+      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 flow to Q:
+      LAPACKE_dorgqr(LAPACK_COL_MAJOR, MATSIZE, MATSIZE, MATSIZE, flow, MATSIZE, tmp11);
+
+      // reset prevtime
+      prevtime=time;
+    }
+
+    // catch abort signal
+    if (g_abort){
+      // print u to stderr if no savefile
+      if (savefile==NULL){
+	savefile=stderr;
+      }
+      break;
+    }
+  }
+
+  if(savefile!=NULL){
+    lyapunov_save_state(flow, u, lyapunov_avg, prevtime, lyapunov_startingtime, savefile, K1, K2, cmd_string, params_string, savefile_string, utfile_string, utfile, COMMAND_LYAPUNOV, algorithm, step, time, nthreads);
+  }
+
+  // save final u to utfile in txt format
+  if(utfile!=NULL){
+    write_vec(u, K1, K2, utfile);
+  }
+
+  lyapunov_free_tmps(lyapunov, lyapunov_avg, flow, u, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9, tmp10, tmp11, fftu1, fftu2, fftu3, fftu4, fft1, ifft, algorithm, algorithm_lyapunov);
+  return(0);
+}
+
+
+// comparison function for qsort
+int compare_double(const void* x, const void* y) {
+  if (*(const double*)x<*(const double*)y) {
+    return(-1);
+  } else if (*(const double*)x>*(const double*)y) {
+    return(+1);
+  } else {
+    return(0);
+  }
+}
+
+// compute tangent flow
+int lyapunov_D_step(
+  double* flow,
+  _Complex double* u,
+  int K1,
+  int K2,
+  int N1,
+  int N2,
+  double nu,
+  double delta,
+  unsigned int algorithm,
+  double L,
+  _Complex double* g,
+  fft_vect fftu1,
+  fft_vect fftu2,
+  fft_vect fftu3,
+  fft_vect fftu4,
+  fft_vect fft1,
+  fft_vect ifft,
+  double* tmp1,
+  double* tmp2,
+  double* tmp3,
+  _Complex double* tmp4,
+  bool irreversible
+){
+  int lx,ly;
+  double alpha;
+  int n;
+
+  // compute fft of u for future use
+  lyapunov_fft_u_T(u,K1,K2,N1,N2,fftu1,fftu2,fftu3,fftu4);
+
+  // compute T
+  lyapunov_T(N1,N2,fftu1,fftu2,fftu3,fftu4,ifft);
+  // save to vect
+  for(lx=0;lx<=K1;lx++){
+    for(ly=(lx>0 ? -K2 : 1);ly<=K2;ly++){
+      tmp4[klookup_sym(lx,ly,K2)]=ifft.fft[klookup(lx,ly,N1,N2)];
+    }
+  }
+  
+  //compute alpha
+  if (irreversible) {
+    alpha=nu;
+  } else {
+    alpha=compute_alpha(u,K1,K2,g,L);
+  }
+
+  // loop over columns
+  for(lx=0;lx<=K1;lx++){
+    for(ly=(lx>0 ? -K2 : 1);ly<=K2;ly++){
+      for(n=0;n<=1;n++){
+	// do not use adaptive step algorithms for the tangent flow: it must be locked to the same times as u
+	if(algorithm==ALGORITHM_RK2){
+	  lyapunov_D_step_rk2(flow+(2*klookup_sym(lx,ly,K2)+n)*MATSIZE, u, K1, K2, N1, N2, alpha, delta, L, g, tmp4, fftu1, fftu2, fftu3, fftu4, fft1, ifft, tmp1, tmp2, irreversible);
+	} else if (algorithm==ALGORITHM_RK4) {
+	  lyapunov_D_step_rk4(flow+(2*klookup_sym(lx,ly,K2)+n)*MATSIZE, u, K1, K2, N1, N2, alpha, delta, L, g, tmp4, fftu1, fftu2, fftu3, fftu4, fft1, ifft, tmp1, tmp2, tmp3, irreversible);
+	} else {
+	  fprintf(stderr,"bug: unknown algorithm for tangent flow: %u, contact ian.jauslin@rutgers.edu\n",algorithm);
+	}
+      }
+    }
+  }
+
+  return(0);
+}
+
+// RK 4 algorithm
+int lyapunov_D_step_rk4(
+  double* flow,
+  _Complex double* u,
+  int K1,
+  int K2,
+  int N1,
+  int N2,
+  double nu,
+  double delta,
+  double L,
+  _Complex double* g,
+  _Complex double* T,
+  fft_vect fftu1,
+  fft_vect fftu2,
+  fft_vect fftu3,
+  fft_vect fftu4,
+  fft_vect fft1,
+  fft_vect ifft,
+  double* tmp1,
+  double* tmp2,
+  double* tmp3,
+  bool irreversible
+){
+  int i;
+
+  // k1
+  lyapunov_D_rhs(tmp1, flow, u, K1, K2, N1, N2, nu, L, g, T, fftu1, fftu2, fftu3, fftu4, fft1, ifft, irreversible);
+  // add to output
+  for(i=0;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 ** 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
+){
+  // 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));
+  *flow=calloc(MATSIZE*MATSIZE,sizeof(double));
+
+  return(0);
+}
+
+// release vectors
+int lyapunov_free_tmps(
+  double * lyapunov,
+  double * lyapunov_avg,
+  double * flow,
+  _Complex double * u,
+  double * tmp1,
+  double * tmp2,
+  double * tmp3,
+  _Complex double * tmp4,
+  _Complex double * tmp5,
+  _Complex double * tmp6,
+  _Complex double * tmp7,
+  _Complex double * tmp8,
+  _Complex double * tmp9,
+  _Complex double * tmp10,
+  double * tmp11,
+  fft_vect fftu1,
+  fft_vect fftu2,
+  fft_vect fftu3,
+  fft_vect fftu4,
+  fft_vect fft1,
+  fft_vect ifft,
+  unsigned int algorithm,
+  unsigned int algorithm_lyapunov
+){
+  free(flow);
+  free(lyapunov);
+  free(lyapunov_avg);
+
+  fftw_destroy_plan(fftu2.fft_plan);
+  fftw_destroy_plan(fftu3.fft_plan);
+  fftw_destroy_plan(fftu4.fft_plan);
+  fftw_free(fftu2.fft);
+  fftw_free(fftu3.fft);
+  fftw_free(fftu4.fft);
+
+  ns_free_tmps(u, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9, tmp10, fft1, fftu1, ifft, algorithm);
+
+  free(tmp11);
+  if(algorithm_lyapunov==ALGORITHM_RK2){
+    free(tmp1);
+    free(tmp2);
+  } else if (algorithm_lyapunov==ALGORITHM_RK4){
+    free(tmp1);
+    free(tmp2);
+    free(tmp3);
+  } else {
+    fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm_lyapunov);
+  };
+
+  return(0);
+}
+
+// save state of lyapunov computation
+int lyapunov_save_state(
+  double* flow,
+  _Complex double* u,
+  double* lyapunov_avg,
+  double prevtime,
+  double lyapunov_startingtime,
+  FILE* savefile,
+  int K1,
+  int K2,
+  const char* cmd_string,
+  const char* params_string,
+  const char* savefile_string,
+  const char* utfile_string,
+  FILE* utfile,
+  unsigned int command,
+  unsigned int algorithm,
+  double step,
+  double time,
+  unsigned int nthreads
+){
+  // save u and step
+  save_state(u, savefile, K1, K2, cmd_string, params_string, savefile_string, utfile_string, utfile, command, algorithm, step, time, nthreads);
+
+  if(savefile!=stderr && savefile!=stdout){
+    // save tangent flow
+    write_mat2_bin(flow,K1,K2,savefile);
+    // save time of QR decomposition
+    fwrite(&prevtime, sizeof(double), 1, savefile);
+    // save time at which the lyapunov computation started
+    fwrite(&lyapunov_startingtime, sizeof(double), 1, savefile);
+    // save lyapunov averages
+    write_vec2_bin(lyapunov_avg,K1,K2,savefile);
+  }
+
+  return 0;
+}
--- a/src/lyapunov.h
+++ b/src/lyapunov.h
@ -0,0 +1,64 @@
+/*
+Copyright 2017-2025 Ian Jauslin
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+#ifndef LYAPUNOV_H
+#define LYAPUNOV_H
+
+#include "navier-stokes.h"
+
+// compute Lyapunov exponents
+int lyapunov( int K1, int K2, int N1, int N2, double final_time, double lyapunov_reset, unsigned int lyapunov_trigger, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_cost, _Complex double* u0, _Complex double* g, bool irreversible, bool keep_en_cst, double target_en, unsigned int algorithm, unsigned int algorithm_lyapunov, double starting_time, unsigned int nthreads, double* flow0, double* lyapunov_avg0, double prevtime, 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);
+
+// 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);
+
+// release vectors
+int lyapunov_free_tmps( double * lyapunov, double * lyapunov_avg, double * flow, _Complex double * u, double * tmp1, double * tmp2, double * tmp3, _Complex double * tmp4, _Complex double * tmp5, _Complex double * tmp6, _Complex double * tmp7, _Complex double * tmp8, _Complex double * tmp9, _Complex double * tmp10, double * tmp11, fft_vect fftu1, fft_vect fftu2, fft_vect fftu3, fft_vect fftu4, fft_vect fft1, fft_vect ifft, unsigned int algorithm, unsigned int algorithm_lyapunov);
+
+// 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-2023 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,21 +16,28 @@ limitations under the License.

 #define VERSION "0.1"

+#include <math.h>
 #include <signal.h>
 #include <string.h>
 #include <stdlib.h>
 #include <errno.h>
+#include <wordexp.h>

 #include "constants.cpp"
 #include "driving.h"
+#include "dstring.h"
 #include "init.h"
 #include "int_tools.h"
+#include "io.h"
+#include "lyapunov.h"
 #include "navier-stokes.h"


 // structure to store parameters, to make it easier to pass parameters to CLI functions
 typedef struct nstrophy_parameters {
-  double init_en; // initial value for the energy for ins and enstrophy for rns
+  double init_energy;
+  double init_enstrophy;
+  unsigned int init_enstrophy_or_energy; //whether to fix the initial enstrophy or energy
  bool irreversible;
  int K1;
  int K2;
@ -43,16 +50,20 @@ typedef struct nstrophy_parameters {
  double adaptive_tolerance;
  double adaptive_factor;
  double max_delta;
-  unsigned int adaptive_norm;
+  unsigned int adaptive_cost;
  double print_freq;
  int seed;
  double starting_time;
  unsigned int driving;
  unsigned int init;
  unsigned int algorithm;
+  unsigned int algorithm_lyapunov;
  bool keep_en_cst;
  FILE* initfile;
  FILE* drivingfile;
+  double lyapunov_reset;
+  unsigned int lyapunov_trigger;
+  bool print_alpha;
 } nstrophy_parameters;

 // usage message
@ -61,14 +72,19 @@ int print_usage();
 int print_params(nstrophy_parameters parameters, char* initfile_str, char* drivingfile_str, FILE* file);

 // read command line arguments
-int read_args(int argc, const char* argv[], char** params, unsigned int* command, unsigned int* nthreads, char** savefile_str, char** utfile_str);
-int read_params(char* param_str, nstrophy_parameters* parameters, char** initfile_str, char** drivingfile_str);
-int set_parameter(char* lhs, char* rhs, nstrophy_parameters* parameters, bool* setN1, bool* setN2, char** initfile_str, char** drivingfile_str);
+int read_args(int argc, const char* argv[], dstring* params, unsigned int* command, unsigned int* nthreads, dstring* savefile_str, dstring* utfile_str, dstring* resumefile_str);
+int set_default_params(nstrophy_parameters* parameters);
+int read_params(dstring param_str, nstrophy_parameters* parameters, dstring* initfile_str, dstring* drivingfile_str);
+int set_parameter(char* lhs, char* rhs, nstrophy_parameters* parameters, bool* setN1, bool* setN2, dstring* initfile_str, dstring* drivingfile_str);
+// read args from file
+int args_from_file(dstring* params, unsigned int* command, unsigned int* nthreads, dstring* savefile_str, dstring* utfile_str, char* file_str);

 // set driving force
 _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);
@ -77,9 +93,17 @@ void sig_handler (int signo);
 volatile bool g_abort = false;
 // signal handler
 void sig_handler (int signo){
-  if (signo == SIGINT || signo == SIGTERM){
+  if (signo == SIGINT){
+    fprintf(stderr,"received signal SIGINT, interrupting computation\n");
    g_abort = true;
  }
+  else if (signo == SIGTERM){
+    fprintf(stderr,"received signal SIGTERM, interrupting computation\n");
+    g_abort = true;
+  }
+  else{
+    fprintf(stderr,"received signal %d\n",signo);
+  }
 }


@ -87,47 +111,121 @@ int main (
  int argc,
  const char* argv[]
 ){
-  char* param_str=NULL;
+  dstring param_str;
  nstrophy_parameters parameters;
  int ret;
  unsigned int command;
  unsigned int nthreads=1;
  _Complex double* u0;
  _Complex double *g;
-  char* savefile_str=NULL;
-  char* utfile_str=NULL;
-  char* initfile_str=NULL;
-  char* drivingfile_str=NULL;
+  double* lyapunov_flow0;
+  double* lyapunov_avg0;
+  double lyapunov_prevtime;
+  double lyapunov_startingtime;
+  dstring savefile_str;
+  dstring utfile_str;
+  dstring initfile_str;
+  dstring drivingfile_str;
+  dstring resumefile_str;
  FILE* savefile=NULL;
  FILE* utfile=NULL;
+  bool resume=false;

  command=0;

+  // init strings
+  dstring_init(&param_str, 64);
+  dstring_init(&savefile_str, 64);
+  dstring_init(&utfile_str, 64);
+  dstring_init(&initfile_str, 64);
+  dstring_init(&drivingfile_str, 64);
+  dstring_init(&resumefile_str, 64);
+
  // read command line arguments
-  ret=read_args(argc, argv, &param_str, &command, &nthreads, &savefile_str, &utfile_str);
+  ret=read_args(argc, argv, &param_str, &command, &nthreads, &savefile_str, &utfile_str, &resumefile_str);
  if(ret<0){
-    return(-1);
+    dstring_free(param_str);
+    dstring_free(savefile_str);
+    dstring_free(utfile_str);
+    dstring_free(initfile_str);
+    dstring_free(drivingfile_str);
+    dstring_free(resumefile_str);
+    return(ret);
  }

+  // set default params
+  set_default_params(&parameters);
  // read params
  ret=read_params(param_str, &parameters, &initfile_str, &drivingfile_str);
  if(ret<0){
-    return(-1);
+    dstring_free(param_str);
+    dstring_free(savefile_str);
+    dstring_free(utfile_str);
+    dstring_free(initfile_str);
+    dstring_free(drivingfile_str);
+    dstring_free(resumefile_str);
+    return(ret);
  }

+  // 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);
+      dstring_free(savefile_str);
+      dstring_free(utfile_str);
+      dstring_free(initfile_str);
+      dstring_free(drivingfile_str);
+      dstring_free(resumefile_str);
+      return(ret);
+    }
+    // read params
+    ret=read_params(param_str, &parameters, &initfile_str, &drivingfile_str);
+    if(ret<0){
+      dstring_free(param_str);
+      dstring_free(savefile_str);
+      dstring_free(utfile_str);
+      dstring_free(initfile_str);
+      dstring_free(drivingfile_str);
+      dstring_free(resumefile_str);
+      return(ret);
+    }
+
+    // reread arguments (to allow overrides from the command line, but do not override the command)
+    unsigned int dummy_command;
+    read_args(argc, argv, &param_str, &dummy_command, &nthreads, &savefile_str, &utfile_str, &resumefile_str);
+    // reread params
+    read_params(param_str, &parameters, &initfile_str, &drivingfile_str);
+  }
+
+  // free strings
+  dstring_free(resumefile_str);
+
  // open initfile
-  if(initfile_str!=NULL){
-    parameters.initfile=fopen(initfile_str,"r");
+  if(initfile_str.length!=0){
+    parameters.initfile=fopen(dstring_to_str_noinit(&initfile_str),"r");
    if(parameters.initfile==NULL){
-      fprintf(stderr,"Error opening file '%s' for reading: %s\n", initfile_str, strerror(errno));
+      fprintf(stderr,"Error opening file '%s' for reading: %s\n", dstring_to_str_noinit(&initfile_str), strerror(errno));
+      dstring_free(param_str);
+      dstring_free(savefile_str);
+      dstring_free(utfile_str);
+      dstring_free(initfile_str);
+      dstring_free(drivingfile_str);
      return(-1);
    }
  }
  // open drivingfile
-  if(drivingfile_str!=NULL){
-    parameters.drivingfile=fopen(drivingfile_str,"r");
+  if(drivingfile_str.length!=0){
+    parameters.drivingfile=fopen(dstring_to_str_noinit(&drivingfile_str),"r");
    if(parameters.drivingfile==NULL){
-      fprintf(stderr,"Error opening file '%s' for reading: %s\n", drivingfile_str, strerror(errno));
+      fprintf(stderr,"Error opening file '%s' for reading: %s\n", dstring_to_str_noinit(&drivingfile_str), strerror(errno));
+      dstring_free(param_str);
+      dstring_free(savefile_str);
+      dstring_free(utfile_str);
+      dstring_free(initfile_str);
+      dstring_free(drivingfile_str);
      return(-1);
    }
  }
@ -136,6 +234,29 @@ 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){
+    parameters.init_enstrophy=compute_enstrophy(u0, parameters.K1, parameters.K2, parameters.L);
+  }
+  // same with init_energy
+  if (parameters.init_enstrophy_or_energy!=FIX_ENERGY){
+    parameters.init_energy=compute_energy(u0, parameters.K1, parameters.K2);
+  }

  // close initfile (do this early, so that it is possible to use the same file for init and save)
  if (parameters.initfile!=NULL){
@ -147,48 +268,71 @@ int main (
  }

  // open savefile (do this after closing init file)
-  if(savefile_str!=NULL){
-    savefile=fopen(savefile_str,"w");
+  if(savefile_str.length!=0){
+    savefile=fopen(dstring_to_str_noinit(&savefile_str),"w");
    if(savefile==NULL){
-      fprintf(stderr,"Error opening file '%s' for writing: %s\n", savefile_str, strerror(errno));
+      fprintf(stderr,"Error opening file '%s' for writing: %s\n", dstring_to_str_noinit(&savefile_str), strerror(errno));
+      dstring_free(param_str);
+      dstring_free(savefile_str);
+      dstring_free(utfile_str);
+      dstring_free(initfile_str);
+      dstring_free(drivingfile_str);
+      free(g);
+      free(u0);
+      if (command==COMMAND_LYAPUNOV){
+	free(lyapunov_flow0);
+	free(lyapunov_avg0);
+      }
      return(-1);
    }
  }

  // open utfile (do this after closing init file)
-  if(utfile_str!=NULL){
-    utfile=fopen(utfile_str,"w");
+  if(utfile_str.length!=0){
+    utfile=fopen(dstring_to_str_noinit(&utfile_str),"w");
    if(utfile==NULL){
-      fprintf(stderr,"Error opening file '%s' for writing: %s\n", utfile_str, strerror(errno));
+      fprintf(stderr,"Error opening file '%s' for writing: %s\n", dstring_to_str_noinit(&utfile_str), strerror(errno));
+      dstring_free(param_str);
+      dstring_free(savefile_str);
+      dstring_free(utfile_str);
+      dstring_free(initfile_str);
+      dstring_free(drivingfile_str);
+      free(g);
+      free(u0);
+      if (command==COMMAND_LYAPUNOV){
+	free(lyapunov_flow0);
+	free(lyapunov_avg0);
+      }
      return(-1);
    }
  }

  // print parameters
-  print_params(parameters, initfile_str, drivingfile_str, stderr);
-  print_params(parameters, initfile_str, drivingfile_str, stdout);
+  print_params(parameters, dstring_to_str_noinit(&initfile_str), dstring_to_str_noinit(&drivingfile_str), stderr);
+  print_params(parameters, dstring_to_str_noinit(&initfile_str), dstring_to_str_noinit(&drivingfile_str), stdout);

-  // free initfile_str
-  if(initfile_str!=NULL){
-    free(initfile_str);
-  }
-  // free drivingfile_str
-  if(drivingfile_str!=NULL){
-    free(drivingfile_str);
-  }
+  // free strings
+  dstring_free(initfile_str);
+  dstring_free(drivingfile_str);

  // 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_norm, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_en, 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
    signal(SIGINT, sig_handler);
    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_norm, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_en, parameters.algorithm, parameters.print_freq, parameters.starting_time, nthreads, savefile, utfile, (char*)argv[0], param_str, savefile_str, utfile_str);
+    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_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_norm, parameters.starting_time, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_en, parameters.algorithm, nthreads, savefile);
+    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){
+    // 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");
@ -197,6 +341,15 @@ int main (

  free(g);
  free(u0);
+  if (command==COMMAND_LYAPUNOV){
+    free(lyapunov_flow0);
+    free(lyapunov_avg0);
+  }
+
+  // free strings
+  dstring_free(param_str);
+  dstring_free(savefile_str);
+  dstring_free(utfile_str);

  // close savefile
  if (savefile!=NULL){
@ -213,7 +366,7 @@ int main (

 // usage message
 int print_usage(){
-  fprintf(stderr, "usage:\n       nstrophy [-t nthreads] [-p parameters] [-s savefile] [-u u_outfile] <command>\n\n");
+  fprintf(stderr, "usage:\n       nstrophy [-t nthreads] [-p parameters] [-s savefile] [-u u_outfile] <command>\n  where <command> is one of\n    enstrophy\n    uk\n    quiet\n    resume <resume_file>\n\n");
  return(0);
 }

@ -232,7 +385,12 @@ int print_params(
    fprintf(file,"equation=reversible");
  }

-  fprintf(file,", K1=%d, K2=%d, N1=%d, N2=%d, nu=%.15e, delta=%.15e, L=%.15e, init_en=%.15e", parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nu, parameters.delta, parameters.L, parameters.init_en);
+  fprintf(file,", K1=%d, K2=%d, N1=%d, N2=%d, nu=%.15e, delta=%.15e, L=%.15e", parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nu, parameters.delta, parameters.L);
+  if (parameters.init_enstrophy_or_energy==FIX_ENSTROPHY){
+    fprintf(file,", init_enstrophy=%.15e", parameters.init_enstrophy);
+  } else if (parameters.init_enstrophy_or_energy==FIX_ENERGY){
+    fprintf(file,", init_energy=%.15e", parameters.init_energy);
+  }

  switch(parameters.driving){
  case DRIVING_TEST:
@ -292,18 +450,21 @@ int print_params(
  }

  if(parameters.algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){
-    switch(parameters.adaptive_norm){
-    case NORM_L1:
-      fprintf(file,", norm=L1");
+    switch(parameters.adaptive_cost){
+    case COST_L1:
+      fprintf(file,", cost=L1");
      break;
-    case NORM_k3:
-      fprintf(file,", norm=k3");
+    case COST_k3:
+      fprintf(file,", cost=k3");
      break;
-    case NORM_k32:
-      fprintf(file,", norm=k32");
+    case COST_k32:
+      fprintf(file,", cost=k32");
      break;
-    case NORM_ENSTROPHY:
-      fprintf(file,", norm=enstrophy");
+    case COST_ENSTROPHY:
+      fprintf(file,", cost=enstrophy");
+      break;
+    case COST_ALPHA:
+      fprintf(file,", cost=alpha");
      break;
    }
  }
@ -319,14 +480,16 @@ int print_params(
 #define CP_FLAG_NTHREADS 2
 #define CP_FLAG_SAVEFILE 3
 #define CP_FLAG_UTFILE 4
+#define CP_FLAG_RESUME 5
 int read_args(
  int argc,
  const char* argv[],
-  char** params,
+  dstring* params,
  unsigned int* command,
  unsigned int* nthreads,
-  char** savefile_str,
-  char** utfile_str
+  dstring* savefile_str,
+  dstring* utfile_str,
+  dstring* resumefile_str
 ){
  int i;
  int ret;
@ -363,7 +526,7 @@ int read_args(
    }
    // params
    else if(flag==CP_FLAG_PARAMS){
-      *params=(char*)argv[i];
+      str_to_dstring_noinit((char*)argv[i], params);
      flag=0;
    }
    // nthreads
@ -377,12 +540,17 @@ int read_args(
    }
    // savefile
    else if(flag==CP_FLAG_SAVEFILE){
-      *savefile_str=(char*)argv[i];
+      str_to_dstring_noinit((char*)argv[i], savefile_str);
      flag=0;
    }
    // savefile
    else if(flag==CP_FLAG_UTFILE){
-      *utfile_str=(char*)argv[i];
+      str_to_dstring_noinit((char*)argv[i], utfile_str);
+      flag=0;
+    }
+    // resume file
+    else if(flag==CP_FLAG_RESUME){
+      str_to_dstring_noinit((char*)argv[i], resumefile_str);
      flag=0;
    }
    // computation to run
@ -396,38 +564,36 @@ int read_args(
      else if(strcmp(argv[i], "quiet")==0){
 	*command=COMMAND_QUIET;
      }
+      else if(strcmp(argv[i], "resume")==0){
+	*command=COMMAND_RESUME;
+	flag=CP_FLAG_RESUME;
+      }
+      else if(strcmp(argv[i], "lyapunov")==0){
+	*command=COMMAND_LYAPUNOV;
+      }
      else{
 	fprintf(stderr, "error: unrecognized command: '%s'\n",argv[i]);
 	return(-1);
      }
-      flag=0;
    }
  }

+  // check that if the command is 'resume', then resumefile has been set
+  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);
+  }
+
  return(0);
 }

-// read parameters string
-int read_params(
-  char* param_str,
-  nstrophy_parameters* parameters,
-  char** initfile_str,
-  char** drivingfile_str
+// set default parameters
+int set_default_params(
+  nstrophy_parameters* parameters
 ){
-  int ret;
-  // pointer in params
-  char* ptr;
-  // buffer and associated pointer
-  char *buffer_lhs, *lhs_ptr;
-  char *buffer_rhs, *rhs_ptr;
-  // whether N was set explicitly
-  bool setN1=false;
-  bool setN2=false;
-  // whether lhs (false is rhs)
-  bool lhs=true;
-
-  // defaults
-  parameters->init_en=1.54511597324389e+02;
+  // no choice
+  parameters->init_enstrophy_or_energy=0; // do not set init_enstrophy or init_energy here: do that after reading init
  parameters->irreversible=true;
  parameters->K1=16;
  parameters->K2=parameters->K1;
@ -439,7 +605,7 @@ int read_params(
  parameters->adaptive_tolerance=1e-11;
  parameters->adaptive_factor=0.9;
  parameters->max_delta=1e-2;
-  parameters->adaptive_norm=NORM_L1;
+  parameters->adaptive_cost=COST_L1;
  parameters->final_time=100000;
  parameters->print_freq=1;
  parameters->starting_time=0;
@ -450,64 +616,73 @@ int read_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;

-  if (param_str!=NULL){
-    // init
-    buffer_lhs=calloc(sizeof(char),strlen(param_str));
-    lhs_ptr=buffer_lhs;
-    *lhs_ptr='\0';
-    buffer_rhs=calloc(sizeof(char),strlen(param_str));
-    rhs_ptr=buffer_rhs;
-    *rhs_ptr='\0';
+  parameters->print_alpha=false;
+  
+  return(0);
+}

-    for(ptr=param_str;*ptr!='\0';ptr++){
-      switch(*ptr){
-      case '=':
-	// reset buffer
-	rhs_ptr=buffer_rhs;
-	*rhs_ptr='\0';
-	lhs=false;
-	break;
-      case ';':
-	//set parameter
-	ret=set_parameter(buffer_lhs, buffer_rhs, parameters, &setN1, &setN2, initfile_str, drivingfile_str);
-	if(ret<0){
-	  return ret;
-	}
-	// reset buffer
-	lhs_ptr=buffer_lhs;
-	*lhs_ptr='\0';
-	lhs=true;
-	break;
-      default:
-	// add to buffer
-	if (lhs){
-	  *lhs_ptr=*ptr;
-	  lhs_ptr++;
-	  *lhs_ptr='\0';
-	}
-	else{
-	  *rhs_ptr=*ptr;
-	  rhs_ptr++;
-	  *rhs_ptr='\0';
-	}
-	break;
-      }
-    }
+// read parameters string
+int read_params(
+  dstring param_str,
+  nstrophy_parameters* parameters,
+  dstring* initfile_str,
+  dstring* drivingfile_str
+){
+  int ret;
+  unsigned int i;
+  // buffer and associated pointer
+  dstring buffer_lhs;
+  dstring buffer_rhs;
+  // whether N was set explicitly
+  bool setN1=false;
+  bool setN2=false;
+  // which buffer to add to
+  dstring* buffer=&buffer_lhs;

-    // set last param
-    if (*param_str!='\0'){
-      ret=set_parameter(buffer_lhs, buffer_rhs, parameters, &setN1, &setN2, initfile_str, drivingfile_str);
+  // init
+  dstring_init(&buffer_lhs, param_str.length);
+  dstring_init(&buffer_rhs, param_str.length);
+
+  for(i=0;i<param_str.length;i++){
+    switch(param_str.string[i]){
+    case '=':
+      // reset buffer
+      buffer_rhs.length=0;
+      buffer=&buffer_rhs;
+      break;
+    case ';':
+      //set parameter
+      ret=set_parameter(dstring_to_str_noinit(&buffer_lhs), dstring_to_str_noinit(&buffer_rhs), parameters, &setN1, &setN2, initfile_str, drivingfile_str);
      if(ret<0){
 	return ret;
      }
+      // reset buffer
+      buffer_lhs.length=0;
+      buffer=&buffer_lhs;
+      break;
+    default:
+      // add to buffer
+      dstring_append(param_str.string[i],buffer);
+      break;
    }
-
-    // free vects
-    free(buffer_lhs);
-    free(buffer_rhs);
  }

+  // set last param
+  if (param_str.length!=0){
+    ret=set_parameter(dstring_to_str_noinit(&buffer_lhs), dstring_to_str_noinit(&buffer_rhs), parameters, &setN1, &setN2, initfile_str, drivingfile_str);
+    if(ret<0){
+      return ret;
+    }
+  }
+
+  // free vects
+  dstring_free(buffer_lhs);
+  dstring_free(buffer_rhs);
+
  // if N not set explicitly, set it to the smallest power of 2 that is >3*K+1 (the fft is faster on vectors whose length is a power of 2)
  if (!setN1){
    parameters->N1=smallest_pow2(3*(parameters->K1));
@ -516,6 +691,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);
 }

@ -527,8 +708,8 @@ int set_parameter(
  nstrophy_parameters* parameters,
  bool* setN1,
  bool* setN2,
-  char** initfile_str,
-  char** drivingfile_str
+  dstring* initfile_str,
+  dstring* drivingfile_str
 ){
  int ret;

@ -544,10 +725,29 @@ int set_parameter(
      return(-1);
    }
  }
-  else if (strcmp(lhs,"init_en")==0){
-    ret=sscanf(rhs,"%lf",&(parameters->init_en));
+  else if (strcmp(lhs,"init_enstrophy")==0){
+    // check that the energy was not also set
+    if (parameters->init_enstrophy_or_energy==FIX_ENERGY){
+      fprintf(stderr, "error: cannot use init_enstrophy if init_energy has been used\n");
+      return(-1);
+    }
+    parameters->init_enstrophy_or_energy=FIX_ENSTROPHY;
+    ret=sscanf(rhs,"%lf",&(parameters->init_enstrophy));
    if(ret!=1){
-      fprintf(stderr, "error: parameter 'init_en' should be a double\n       got '%s'\n",rhs);
+      fprintf(stderr, "error: parameter 'init_enstrophy' should be a double\n       got '%s'\n",rhs);
+      return(-1);
+    }
+  }
+  else if (strcmp(lhs,"init_energy")==0){
+    // check that the enstrophy was not also set
+    if (parameters->init_enstrophy_or_energy==FIX_ENSTROPHY){
+      fprintf(stderr, "error: cannot use init_energy if init_enstrophy has been used\n");
+      return(-1);
+    }
+    parameters->init_enstrophy_or_energy=FIX_ENERGY;
+    ret=sscanf(rhs,"%lf",&(parameters->init_energy));
+    if(ret!=1){
+      fprintf(stderr, "error: parameter 'init_energy' should be a double\n       got '%s'\n",rhs);
      return(-1);
    }
  }
@ -648,21 +848,24 @@ int set_parameter(
      return(-1);
    }
  }
-  else if (strcmp(lhs,"adaptive_norm")==0){
+  else if (strcmp(lhs,"adaptive_cost")==0){
    if (strcmp(rhs,"L1")==0){
-      parameters->adaptive_norm=NORM_L1;
+      parameters->adaptive_cost=COST_L1;
    }
    else if (strcmp(rhs,"k3")==0){
-      parameters->adaptive_norm=NORM_k3;
+      parameters->adaptive_cost=COST_k3;
    }
    else if (strcmp(rhs,"k32")==0){
-      parameters->adaptive_norm=NORM_k32;
+      parameters->adaptive_cost=COST_k32;
    }
    else if (strcmp(rhs,"enstrophy")==0){
-      parameters->adaptive_norm=NORM_ENSTROPHY;
+      parameters->adaptive_cost=COST_ENSTROPHY;
+    }
+    else if (strcmp(rhs,"alpha")==0){
+      parameters->adaptive_cost=COST_ALPHA;
    }
    else{
-      fprintf(stderr, "error: unrecognized adaptive_norm '%s'\n",rhs);
+      fprintf(stderr, "error: unrecognized adaptive_cost '%s'\n",rhs);
      return(-1);
    }
  }
@ -697,14 +900,12 @@ int set_parameter(
    // matches any argument that starts with 'file:'
    else if (strncmp(rhs,"file:",5)==0){
      parameters->driving=DRIVING_FILE;
-      *drivingfile_str=calloc(sizeof(char), strlen(rhs)-5+1);
-      strcpy(*drivingfile_str, rhs+5);
+      str_to_dstring_noinit(rhs+5,drivingfile_str);
    }
    // matches any argument that starts with 'file_txt:'
    else if (strncmp(rhs,"file_txt:",9)==0){
      parameters->driving=DRIVING_FILE_TXT;
-      *drivingfile_str=calloc(sizeof(char), strlen(rhs)-9+1);
-      strcpy(*drivingfile_str, rhs+9);
+      str_to_dstring_noinit(rhs+9,drivingfile_str);
    }
    else{
      fprintf(stderr, "error: unrecognized driving force '%s'\n",rhs);
@ -722,14 +923,12 @@ int set_parameter(
    // matches any argument that starts with 'file:'
    else if (strncmp(rhs,"file:",5)==0){
      parameters->init=INIT_FILE;
-      *initfile_str=calloc(sizeof(char), strlen(rhs)-5+1);
-      strcpy(*initfile_str, rhs+5);
+      str_to_dstring_noinit(rhs+5,initfile_str);
    }
    // matches any argument that starts with 'file_txt:'
    else if (strncmp(rhs,"file_txt:",9)==0){
      parameters->init=INIT_FILE_TXT;
-      *initfile_str=calloc(sizeof(char), strlen(rhs)-9+1);
-      strcpy(*initfile_str, rhs+9);
+      str_to_dstring_noinit(rhs+9,initfile_str);
    }
    else{
      fprintf(stderr, "error: unrecognized initial condition '%s'\n",rhs);
@ -767,6 +966,56 @@ int set_parameter(
    }
    parameters->keep_en_cst=(tmp==1);
  }
+  else if (strcmp(lhs,"print_alpha")==0){
+    int tmp;
+    ret=sscanf(rhs,"%d",&tmp);
+    if(ret!=1 || (tmp!=0 && tmp!=1)){
+      fprintf(stderr, "error: parameter 'print_alpha' should be 0 or 1\n       got '%s'\n",rhs);
+      return(-1);
+    }
+    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);
@ -775,11 +1024,72 @@ int set_parameter(
  return(0);
 }

+// read args from file
+int args_from_file(
+  dstring* params,
+  unsigned int* command,
+  unsigned int* nthreads,
+  dstring* savefile_str,
+  dstring* utfile_str,
+  char* file_str
+){
+  dstring line;
+  char c;
+
+  // open file
+  FILE* file=fopen(file_str,"r");
+  if(file==NULL){
+    fprintf(stderr,"Error opening file '%s' for reading: %s\n", file_str, strerror(errno));
+    return(-1);
+  }
+
+  // allocate line buffer
+  dstring_init(&line,64);
+
+  while(1){
+    c=fgetc(file);
+
+    // end of file
+    if (feof(file)){
+      break;
+    }
+
+    // newline: read line and reset buffer
+    if(c=='\n' || c=='\r'){
+      // read entry
+      // ignore lines with fewer than three characters
+      if(line.length>=3){
+	// find line starting with "#!"
+	if(line.string[0]=='#' && line.string[1]=='!'){
+	  wordexp_t pwordexp;
+	  wordexp(line.string+2,&pwordexp,0);
+	  // read arguments
+	  read_args(pwordexp.we_wordc, (const char **)pwordexp.we_wordv, params, command, nthreads, savefile_str, utfile_str, NULL);
+	  wordfree(&pwordexp);
+	  // exit
+	  break;
+	}
+      }
+      // reset buffer
+      line.length=0;
+    }
+    // add to buffer
+    else {
+      dstring_append(c,&line);
+    }
+  }
+  dstring_free(line);
+  // close file
+  fclose(file);
+
+  return(0);
+}
+
 // set driving force
 _Complex double* set_driving(
  nstrophy_parameters parameters
 ){
-  _Complex double* g=calloc(sizeof(_Complex double),parameters.K1*(2*parameters.K2+1)+parameters.K2);
+  _Complex double* g=calloc(parameters.K1*(2*parameters.K2+1)+parameters.K2,sizeof(_Complex double));

  switch(parameters.driving){
  case DRIVING_ZERO:
@ -808,25 +1118,19 @@ _Complex double* set_driving(
 _Complex double* set_init(
  nstrophy_parameters* parameters
 ){
-  _Complex double* u0=calloc(sizeof(_Complex double),parameters->K1*(2*parameters->K2+1)+parameters->K2);
+  _Complex double* u0=calloc(parameters->K1*(2*parameters->K2+1)+parameters->K2,sizeof(_Complex double));

  switch(parameters->init){
  case INIT_RANDOM:
-    init_random(u0, parameters->init_en, parameters->K1, parameters->K2, parameters->L, parameters->seed, parameters->irreversible);
+    init_random(u0, parameters->K1, parameters->K2, parameters->seed);
    break;

  case INIT_GAUSSIAN:
-    init_gaussian(u0, parameters->init_en, parameters->K1, parameters->K2, parameters->L, parameters->irreversible);
+    init_gaussian(u0, parameters->K1, parameters->K2);
    break;

  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:
@ -834,9 +1138,76 @@ _Complex double* set_init(
    break;

  default:
-    init_gaussian(u0, parameters->init_en, parameters->K1, parameters->K2, parameters->L, parameters->irreversible);
+    init_gaussian(u0, parameters->K1, parameters->K2);
    break;
  }
+
+  // rescale energy or enstrophy
+  if (parameters->init_enstrophy_or_energy==FIX_ENERGY) {
+    // fix the energy
+    double rescale=compute_energy(u0, parameters->K1, parameters->K2);
+    int kx,ky;
+    for(kx=0;kx<=parameters->K1;kx++){
+      for(ky=(kx>0 ? -parameters->K2 : 1);ky<=parameters->K2;ky++){
+	u0[klookup_sym(kx,ky,parameters->K2)]=u0[klookup_sym(kx,ky,parameters->K2)]*sqrt(parameters->init_energy/rescale);
+      }
+    }
+  } else if (parameters->init_enstrophy_or_energy==FIX_ENSTROPHY) {
+    // fix the enstrophy
+    double rescale=compute_enstrophy(u0, parameters->K1, parameters->K2, parameters->L);
+    int kx,ky;
+    for(kx=0;kx<=parameters->K1;kx++){
+      for(ky=(kx>0 ? -parameters->K2 : 1);ky<=parameters->K2;ky++){
+	u0[klookup_sym(kx,ky,parameters->K2)]=u0[klookup_sym(kx,ky,parameters->K2)]*sqrt(parameters->init_enstrophy/rescale);
+      }
+    }
+  }
  
  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
--- a/src/navier-stokes.h
+++ b/src/navier-stokes.h
@ -1,5 +1,5 @@
 /*
-Copyright 2017-2023 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,8 +21,6 @@ limitations under the License.
 #include <stdbool.h>
 #include <fftw3.h>

-#define M_PI 3.14159265358979323846
-
 // abort signal
 extern volatile bool g_abort;

@ -33,13 +31,13 @@ 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_norm, _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_norm, _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 nthreads, FILE* savefile, FILE* utfile, char* cmd_string, char* params_string, char* savefile_string, char* utfile_string);
+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 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_norm, 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);
+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);


 // initialize vectors for computation
@ -51,16 +49,21 @@ int ns_free_tmps( _Complex double* u, _Complex double* tmp1, _Complex double *tm
 int copy_u( _Complex double* u, _Complex double* u0, int K1, int K2);

 // next time step for Irreversible/reversible Navier-Stokes equation
+int ns_step( unsigned int algorithm, _Complex double* u, int K1, int K2, int N1, int N2, double nu, double* delta, double* next_delta, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_cost, double L, _Complex double* g, double time, double starting_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, bool irreversible, bool keep_en_cst, double target_en);
 // RK4
 int ns_step_rk4( _Complex double* u, int K1, int K2, int N1, int N2, double nu, double delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2,fft_vect ifft, _Complex double* tmp1, _Complex double *tmp2, _Complex double *tmp3, bool irreversible, bool keep_en_cst, double target_en);
 // RK2
 int ns_step_rk2( _Complex double* u, int K1, int K2, int N1, int N2, double nu, double delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2,fft_vect ifft, _Complex double* tmp1, _Complex double *tmp2, bool irreversible, bool keep_en_cst, double target_en);
 // Runge-Kutta-Fehlberg
-int ns_step_rkf45( _Complex double* u, double tolerance, double factor, double max_delta, unsigned int adaptive_norm, int K1, int K2, int N1, int N2, double nu, double* delta, double* next_delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double* k1, _Complex double* k2, _Complex double* k3, _Complex double* k4, _Complex double* k5, _Complex double* k6, _Complex double* tmp, bool irreversible, bool keep_en_cst, double target_en, bool compute_k1);
+int ns_step_rkf45( _Complex double* u, double tolerance, double factor, double max_delta, unsigned int adaptive_cost, int K1, int K2, int N1, int N2, double nu, double* delta, double* next_delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double* k1, _Complex double* k2, _Complex double* k3, _Complex double* k4, _Complex double* k5, _Complex double* k6, _Complex double* tmp, bool irreversible, bool keep_en_cst, double target_en, bool compute_k1);
 // Runge-Kutta-Dromand-Prince
-int ns_step_rkdp54( _Complex double* u, double tolerance, double factor, double max_delta, unsigned int adaptive_norm, int K1, int K2, int N1, int N2, double nu, double* delta, double* next_delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double** k1, _Complex double** k2, _Complex double* k3, _Complex double* k4, _Complex double* k5, _Complex double* k6, _Complex double* tmp, bool irreversible, bool keep_en_cst, double target_en, bool compute_k1);
+int ns_step_rkdp54( _Complex double* u, double tolerance, double factor, double max_delta, unsigned int adaptive_cost, int K1, int K2, int N1, int N2, double nu, double* delta, double* next_delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double** k1, _Complex double** k2, _Complex double* k3, _Complex double* k4, _Complex double* k5, _Complex double* k6, _Complex double* tmp, bool irreversible, bool keep_en_cst, double target_en, bool compute_k1);
 // Runge-Kutta-Bogacki-Shampine
-int ns_step_rkbs32( _Complex double* u, double tolerance, double factor, double max_delta, unsigned int adaptive_norm, int K1, int K2, int N1, int N2, double nu, double* delta, double* next_delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double** k1, _Complex double* k2, _Complex double* k3, _Complex double** k4, _Complex double* tmp, bool irreversible, bool keep_en_cst, double target_en, bool compute_k1);
+int ns_step_rkbs32( _Complex double* u, double tolerance, double factor, double max_delta, unsigned int adaptive_cost, int K1, int K2, int N1, int N2, double nu, double* delta, double* next_delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double** k1, _Complex double* k2, _Complex double* k3, _Complex double** k4, _Complex double* tmp, bool irreversible, bool keep_en_cst, double target_en, bool compute_k1);
+
+// cost function for the adaptive iterations
+double ns_adaptive_cost( _Complex double* u, _Complex double* U, unsigned int adaptive_cost, int K1, int K2, _Complex double* g, double L);
+

 // right side of Irreversible/reversible Navier-Stokes equation
 int ns_rhs( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2, double nu, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft, bool irreversible);
@ -71,6 +74,13 @@ int ns_T( _Complex double* u, int K1, int K2, int N1, int N2, fft_vect fft1, fft
 // convolution term in right side of equation (computed without fft)
 int ns_T_nofft( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2);

+/*
+// Jacobian of rhs
+int ns_D_rhs( _Complex double* out, _Complex double* u, int K1, int K2, double nu, double L, bool irreversible);
+// derivative of T with respect to u_k
+_Complex double ns_d_T( _Complex double* u, int kx, int ky, int lx, int ly, int K2);
+*/
+
 // compute alpha
 double compute_alpha( _Complex double* u, int K1, int K2, _Complex double* g, double L);
 // compute energy
Author	SHA1	Message	Date
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
Ian Jauslin	9ec233fed5	Replace init_en with init_enstrophy and init_energy	2024-12-11 15:48:11 -05:00
Ian Jauslin	d81ad18618	Rewrite cost function for adaptive step	2024-11-18 19:11:00 -05:00
Ian Jauslin	9fa10c8db4	print_alpha	2024-11-07 14:39:03 -05:00
Ian Jauslin	cf48b23d4d	Update copyright	2024-10-15 11:47:13 -04:00
Ian Jauslin	c7ff1384e1	Update README	2024-10-15 11:45:28 -04:00
Ian Jauslin	0404bd0cf0	Print negative values of alpha	2024-10-15 11:41:31 -04:00
Ian Jauslin	7675148e7d	Change order of outputs in enstrophy	2024-10-15 11:36:48 -04:00
Ian Jauslin	f93ada07f0	Fix order of arguments in calloc	2024-10-15 11:31:53 -04:00
Ian Jauslin	41ff1919f0	Change outputs of enstrophy	2024-10-15 11:28:07 -04:00
Ian Jauslin	41a5a4ba3f	save_state to its own function	2024-03-13 11:24:47 +01:00
Ian Jauslin	0599f69dc7	Fix lyapunov computation	2024-03-12 01:53:48 +01:00
Ian Jauslin	198867751d	Split real and imag part in jacobian	2024-03-11 23:54:11 +01:00
Ian Jauslin	9059a24c23	Fix MATSIZE	2024-02-20 10:55:15 -05:00
Ian Jauslin	f9171aa355	Fix average	2024-02-20 10:46:53 -05:00
Ian Jauslin	c284587105	Sort Lyapunov exponents and print more	2024-02-20 09:32:32 -05:00
Ian Jauslin	bf8e4b728e	Fix lyapunov	2024-02-19 22:05:30 -05:00
Ian Jauslin	ec98098709	Lyapunov parameters	2024-02-19 21:39:48 -05:00
Ian Jauslin	1c3f4b9f2f	Print Lyapunov exponents	2024-02-19 19:09:26 -05:00
Ian Jauslin	15bcb8b467	Lyapunov command in main	2024-02-19 19:05:35 -05:00
Ian Jauslin	b31fab8eae	Computation of Lyapunov exponents	2024-02-19 19:01:31 -05:00
Ian Jauslin	89601d19d1	Jacobian of un->un+1	2024-02-19 15:20:19 -05:00
Ian Jauslin	9ecf4413a5	write ns_step routine that handles all the algorithms	2024-02-19 14:43:04 -05:00
Ian Jauslin	9fe9e3d96f	Jacobian of rhs	2024-02-19 12:20:21 -05:00
Ian Jauslin	63e983b460	Link openblas	2024-02-19 11:41:47 -05:00
Ian Jauslin	97a127a8db	Document Lyapunov exponents	2024-02-19 11:41:35 -05:00
Ian Jauslin	9b9734c4c2	switch point to subsubsection	2024-02-18 14:18:02 -05:00
Ian Jauslin	6c12e47105	Do not override command when using 'resume'	2023-11-03 18:02:12 -04:00
Ian Jauslin	209ba06cbf	memory fix	2023-11-03 17:41:30 -04:00
Ian Jauslin	16c80d2305	strlen in dstrings.c	2023-11-03 16:48:18 -04:00
Ian Jauslin	f9aac70796	dstring in read_args_from_file	2023-11-03 16:45:51 -04:00
Ian Jauslin	9d232a8fca	dstrings in read_params	2023-11-03 16:42:27 -04:00
Ian Jauslin	04a15dd2c7	Replace argv strings with dstrings	2023-11-03 16:04:53 -04:00
Ian Jauslin	dd9bd74c83	dynamic strings	2023-11-03 12:25:02 -04:00
Ian Jauslin	3d94694017	Resume command Not functional yet: need to deal with memory of argv	2023-11-03 12:08:25 -04:00