From b31fab8eaed204c87c5e14c037434e4a5b996727 Mon Sep 17 00:00:00 2001
From: Ian Jauslin <ian.jauslin@rutgers.edu>
Date: Mon, 19 Feb 2024 19:01:31 -0500
Subject: [PATCH] Computation of Lyapunov exponents

---
 src/lyapunov.c | 208 ++++++++++++++++++++++++++++++++++++++++++++++++-
 src/lyapunov.h |  10 +++
 2 files changed, 215 insertions(+), 3 deletions(-)

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