From 89601d19d1437aa5f4b711ca1a4ea86fdd581435 Mon Sep 17 00:00:00 2001
From: Ian Jauslin <ian.jauslin@rutgers.edu>
Date: Mon, 19 Feb 2024 15:20:19 -0500
Subject: [PATCH] Jacobian of un->un+1

---
 src/constants.cpp   |   1 +
 src/lyapunov.c      | 107 ++++++++++++++++++++++++++++++++++++++++++++
 src/lyapunov.h      |  25 +++++++++++
 src/navier-stokes.c |   1 +
 src/navier-stokes.h |   2 -
 5 files changed, 134 insertions(+), 2 deletions(-)
 create mode 100644 src/lyapunov.c
 create mode 100644 src/lyapunov.h

diff --git a/src/constants.cpp b/src/constants.cpp
index 21620c9..8a0f054 100644
--- a/src/constants.cpp
+++ b/src/constants.cpp
@@ -14,6 +14,7 @@ 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
diff --git a/src/lyapunov.c b/src/lyapunov.c
new file mode 100644
index 0000000..28ded4f
--- /dev/null
+++ b/src/lyapunov.c
@@ -0,0 +1,107 @@
+/*
+Copyright 2017-2023 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 <cblas.h>
+
+#define MATSIZE (K1*(2*(K2+1)+K2))
+
+// Jacobian of u_n -> u_{n+1}
+int ns_D_step(
+  _Complex double* out,
+  _Complex double* un,
+  _Complex double* un_next,
+  int K1,
+  int K2,
+  int N1,
+  int N2,
+  double nu,
+  double epsilon,
+  double delta,
+  unsigned int algorithm,
+  double adaptive_tolerance,
+  double adaptive_factor,
+  double max_delta,
+  unsigned int adaptive_norm,
+  double L,
+  _Complex double* g,
+  double time,
+  fft_vect fft1,
+  fft_vect fft2,
+  fft_vect ifft,
+  _Complex double* tmp1,
+  _Complex double* tmp2,
+  _Complex double* tmp3,
+  _Complex double* tmp4,
+  _Complex double* tmp5,
+  _Complex double* tmp6,
+  _Complex double* tmp7,
+  _Complex double* tmp8,
+  bool irreversible,
+  bool keep_en_cst,
+  double target_en
+){
+  int lx,ly;
+  int i;
+  double step, next_step;
+
+  for(lx=0;lx<=K1;lx++){
+    for(ly=(lx>0 ? -K2 : 1);ly<=K2;ly++){
+      // derivative in the real direction
+      // finite difference vector
+      for (i=0;i<MATSIZE;i++){
+	if(i==klookup_sym(lx,ly,K2)){
+	  tmp1[i]=un[i]+epsilon;
+	}else{
+	  tmp1[i]=un[i];
+	}
+      }
+      // compute step 
+      // reinitialize step
+      step=delta;
+      next_step=delta;
+      ns_step(algorithm, tmp1, K1, K2, N1, N2, nu, &step, &next_step, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, L, g, time, time, fft1, fft2, ifft, &tmp2, &tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, irreversible, keep_en_cst, target_en);
+      // compute derivative
+      for (i=0;i<MATSIZE;i++){
+	out[i]=(tmp1[i]-un_next[i])/epsilon;
+      }
+
+      // derivative in the imaginary direction
+      // finite difference vector
+      for (i=0;i<MATSIZE;i++){
+	if(i==klookup_sym(lx,ly,K2)){
+	  tmp1[i]=un[i]+epsilon*I;
+	}else{
+	  tmp1[i]=un[i];
+	}
+      }
+      // compute step 
+      // reinitialize step
+      step=delta;
+      next_step=delta;
+      ns_step(algorithm, tmp1, K1, K2, N1, N2, nu, &step, &next_step, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, L, g, time, time, fft1, fft2, ifft, &tmp2, &tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, irreversible, keep_en_cst, target_en);
+      // compute derivative
+      for (i=0;i<MATSIZE;i++){
+	out[i]+=-(tmp1[i]-un_next[i])/epsilon*I;
+      }
+    }
+  }
+
+  return(0);
+}
+
+
diff --git a/src/lyapunov.h b/src/lyapunov.h
new file mode 100644
index 0000000..d36824a
--- /dev/null
+++ b/src/lyapunov.h
@@ -0,0 +1,25 @@
+/*
+Copyright 2017-2023 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"
+
+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);
+#endif
+
+
diff --git a/src/navier-stokes.c b/src/navier-stokes.c
index ef03bf3..28f1fcd 100644
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -1355,6 +1355,7 @@ int ns_T_nofft(
   return 0;
 }
 
+
 /*
 // Jacobian of rhs
 int ns_D_rhs(
diff --git a/src/navier-stokes.h b/src/navier-stokes.h
index ae35bde..68c7caa 100644
--- a/src/navier-stokes.h
+++ b/src/navier-stokes.h
@@ -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;