Compare commits
14 Commits
6c12e47105
...
9059a24c23
Author | SHA1 | Date |
---|---|---|
Ian Jauslin | 9059a24c23 | |
Ian Jauslin | f9171aa355 | |
Ian Jauslin | c284587105 | |
Ian Jauslin | bf8e4b728e | |
Ian Jauslin | ec98098709 | |
Ian Jauslin | 1c3f4b9f2f | |
Ian Jauslin | 15bcb8b467 | |
Ian Jauslin | b31fab8eae | |
Ian Jauslin | 89601d19d1 | |
Ian Jauslin | 9ecf4413a5 | |
Ian Jauslin | 9fe9e3d96f | |
Ian Jauslin | 63e983b460 | |
Ian Jauslin | 97a127a8db | |
Ian Jauslin | 9b9734c4c2 |
2
Makefile
2
Makefile
|
@ -30,7 +30,7 @@ LD=$(CC)
|
|||
#INCLUDES =
|
||||
|
||||
#LIBDIRS =
|
||||
LIBS = -lm -lfftw3 -lfftw3_threads
|
||||
LIBS = -lm -lfftw3 -lfftw3_threads -lopenblas_64
|
||||
|
||||
|
||||
override LDFLAGS +=$(LIBDIRS)$(LIBS)
|
||||
|
|
|
@ -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
|
||||
|
@ -101,9 +103,8 @@ We truncate the Fourier modes and assume that $\hat u_k=0$ if $|k_1|>K_1$ or $|k
|
|||
\mathcal K:=\{(k_1,k_2),\ |k_1|\leqslant K_1,\ |k_2|\leqslant K_2\}
|
||||
.
|
||||
\end{equation}
|
||||
\bigskip
|
||||
|
||||
\point{\bf Runge-Kutta methods}.
|
||||
\subsubsection{Runge-Kutta methods}.
|
||||
To solve the equation numerically, we will use Runge-Kutta methods, which compute an approximate value $\hat u_k^{(n)}$ for $\hat u_k(t_n)$.
|
||||
{\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:
|
||||
|
@ -198,9 +199,8 @@ It can be made by specifying the parameter {\tt adaptive\_norm}.
|
|||
\end{equation}
|
||||
This norm is selected by choosing {\tt adaptive\_norm=enstrophy}.
|
||||
\end{itemize}
|
||||
\bigskip
|
||||
|
||||
\point{\bf Reality}.
|
||||
\subsubsection{Reality}.
|
||||
Since $U$ is real, $\hat U_{-k}=\hat U_k^*$, and so
|
||||
\begin{equation}
|
||||
\hat u_{-k}=\hat u_k^*
|
||||
|
@ -222,9 +222,8 @@ 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$.
|
||||
\bigskip
|
||||
|
||||
\point{\bf FFT}. We compute T using a fast Fourier transform, defined as
|
||||
\subsubsection{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}
|
||||
|
@ -267,9 +266,8 @@ Therefore,
|
|||
\mathcal F\left(q_x|q|\hat u_q\right)(n)
|
||||
\right)(k)
|
||||
\end{equation}
|
||||
\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,18 +368,83 @@ 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}
|
||||
|
||||
\subsubsection{Lyapunov exponents}
|
||||
\indent
|
||||
To compute the Lyapunov exponents, we must first compute the Jacobian of $\hat u^{(n)}\mapsto\hat u^{(n+1)}$.
|
||||
This map is always of Runge-Kutta type, that is,
|
||||
\begin{equation}
|
||||
\hat u^{(n+1)}=\hat u^{(n)}+\delta\sum_{i=1}^q w_i\mathfrak F(\hat u^{(n)})
|
||||
\end{equation}
|
||||
(see\-~(\ref{ins_k})), so
|
||||
\begin{equation}
|
||||
D\hat u^{(n+1)}=\mathds 1+\delta\sum_{i=1}^q w_iD\mathfrak F(\hat u^{(n)})
|
||||
.
|
||||
\end{equation}
|
||||
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
|
||||
|
||||
\point{\bf Numerical instability}.
|
||||
\indent
|
||||
The Lyapunov exponents are then defined as
|
||||
\begin{equation}
|
||||
\frac1n\log\mathrm{spec}(\pi_i)
|
||||
,\quad
|
||||
\pi_i:=\prod_{i=1}^nD\hat u^{(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\in\mathbb N_*$ (set by adjusting the {\tt lyanpunov\_reset} parameter), and, when $n$ is a multiple of $\mathfrak L_r$, we rescale the eigenvalues of $\pi_i$ to 1.
|
||||
To do so, we perform a $QR$ decomposition:
|
||||
\begin{equation}
|
||||
\pi_{\alpha\mathfrak L_r}
|
||||
=Q^{(\alpha)}R^{(\alpha)}
|
||||
\end{equation}
|
||||
where $Q^{(\alpha)}$ is diagonal and $R^{(\alpha)}$ is an orthogonal matrix, and we divide by $Q^{(\alpha)}$ (thus only keeping $R^{(\alpha)}$).
|
||||
Thus, we replace
|
||||
\begin{equation}
|
||||
\pi_{\alpha\mathfrak L_r+i}\mapsto R^{(\alpha)}\prod_{j=1}^iD\hat u^{(j)}
|
||||
.
|
||||
\end{equation}
|
||||
The Lyapunov exponents at time $\alpha\mathfrak L_r$ are then
|
||||
\begin{equation}
|
||||
\frac1{\alpha\mathfrak L_r}\sum_{\beta=1}^\alpha\log\mathrm{spec}(Q^{(\beta)})
|
||||
.
|
||||
\end{equation}
|
||||
|
||||
\subsubsection{Numerical instability}.
|
||||
In order to prevent the algorithm from blowing up, it is necessary to impose the reality of $u(x)$ by hand, otherwise, truncation errors build up, and lead to divergences.
|
||||
It is sufficient to ensure that the convolution term $T(\hat u,k)$ satisfies $T(\hat u,-k)=T(\hat u,k)^*$.
|
||||
After imposing this condition, the algorithm no longer blows up, but it is still unstable (for instance, increasing $K_1$ or $K_2$ leads to very different results).
|
||||
|
|
|
@ -14,11 +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
|
||||
|
|
|
@ -0,0 +1,340 @@
|
|||
/*
|
||||
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 <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[i]=0;
|
||||
}
|
||||
|
||||
// save times at which Lyapunov exponents are computed
|
||||
double prevtime=starting_time;
|
||||
|
||||
// iterate
|
||||
time=starting_time;
|
||||
while(final_time<0 || time<final_time){
|
||||
// copy before step
|
||||
copy_u(prevu, u, K1, K2);
|
||||
prev_step=step;
|
||||
|
||||
ns_step(algorithm, u, K1, K2, N1, N2, nu, &step, &next_step, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, L, g, time, starting_time, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en);
|
||||
|
||||
// compute Jacobian
|
||||
// do not touch tmp1, it might be used in the next step
|
||||
ns_D_step(Du, prevu, u, K1, K2, N1, N2, nu, D_epsilon, prev_step, algorithm, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, L, g, time, fft1, fft2, ifft, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9, irreversible, keep_en_cst, target_en);
|
||||
|
||||
// multiply Jacobian
|
||||
// switch to column major order, so transpose Du
|
||||
cblas_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);
|
||||
}
|
||||
|
||||
// sort lyapunov exponents
|
||||
qsort(lyapunov, MATSIZE, sizeof(double), compare_double);
|
||||
|
||||
// average lyapunov
|
||||
for(i=0; i<MATSIZE; i++){
|
||||
// exclude inf
|
||||
if((! isinf(lyapunov[i])) && (time>starting_time)){
|
||||
lyapunov_avg[i]=lyapunov_avg[i]*(prevtime-starting_time)/(time-starting_time)+lyapunov[i]*(time-prevtime)/(time-starting_time);
|
||||
}
|
||||
}
|
||||
|
||||
// print
|
||||
for(i=0; i<MATSIZE; i++){
|
||||
printf("% .15e % .15e % .15e\n",time, lyapunov[i], lyapunov_avg[i]);
|
||||
}
|
||||
printf("\n");
|
||||
fprintf(stderr,"% .15e",time);
|
||||
// print largest and smallest lyapunov exponent
|
||||
if(MATSIZE>0){
|
||||
fprintf(stderr," % .15e % .15e\n", lyapunov[0], lyapunov[MATSIZE-1]);
|
||||
}
|
||||
|
||||
// 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);
|
||||
}
|
||||
|
||||
// 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);
|
||||
}
|
||||
}
|
||||
|
||||
// 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++){
|
||||
// use row major order
|
||||
out[klookup_sym(lx,ly,K2)*MATSIZE+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++){
|
||||
// use row major order
|
||||
out[klookup_sym(lx,ly,K2)*MATSIZE+i]=-(tmp1[i]-un_next[i])/epsilon*I;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return(0);
|
||||
}
|
||||
|
||||
|
||||
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);
|
||||
}
|
|
@ -0,0 +1,38 @@
|
|||
/*
|
||||
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"
|
||||
|
||||
// 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);
|
||||
|
||||
// comparison function for qsort
|
||||
int compare_double(const void* x, const void* y);
|
||||
|
||||
// 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
|
||||
|
||||
|
23
src/main.c
23
src/main.c
|
@ -27,6 +27,7 @@ limitations under the License.
|
|||
#include "dstring.h"
|
||||
#include "init.h"
|
||||
#include "int_tools.h"
|
||||
#include "lyapunov.h"
|
||||
#include "navier-stokes.h"
|
||||
|
||||
|
||||
|
@ -55,6 +56,8 @@ typedef struct nstrophy_parameters {
|
|||
bool keep_en_cst;
|
||||
FILE* initfile;
|
||||
FILE* drivingfile;
|
||||
double lyapunov_reset;
|
||||
double D_epsilon;
|
||||
} nstrophy_parameters;
|
||||
|
||||
// usage message
|
||||
|
@ -279,6 +282,9 @@ int main (
|
|||
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);
|
||||
}
|
||||
else if(command==COMMAND_LYAPUNOV){
|
||||
lyapunov(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.lyapunov_reset, parameters.nu, parameters.D_epsilon, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_norm, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_en, parameters.algorithm, parameters.starting_time, nthreads);
|
||||
}
|
||||
else if(command==0){
|
||||
fprintf(stderr, "error: no command specified\n");
|
||||
print_usage();
|
||||
|
@ -501,6 +507,9 @@ int read_args(
|
|||
*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);
|
||||
|
@ -786,6 +795,20 @@ int set_parameter(
|
|||
return(-1);
|
||||
}
|
||||
}
|
||||
else if (strcmp(lhs,"lyapunov_reset")==0){
|
||||
ret=sscanf(rhs,"%lf",&(parameters->lyapunov_reset));
|
||||
if(ret!=1){
|
||||
fprintf(stderr, "error: parameter 'lyapunov_reset' should be a double\n got '%s'\n",rhs);
|
||||
return(-1);
|
||||
}
|
||||
}
|
||||
else if (strcmp(lhs,"D_epsilon")==0){
|
||||
ret=sscanf(rhs,"%lf",&(parameters->D_epsilon));
|
||||
if(ret!=1){
|
||||
fprintf(stderr, "error: parameter 'D_epsilon' should be a double\n got '%s'\n",rhs);
|
||||
return(-1);
|
||||
}
|
||||
}
|
||||
else if (strcmp(lhs,"driving")==0){
|
||||
if (strcmp(rhs,"zero")==0){
|
||||
parameters->driving=DRIVING_ZERO;
|
||||
|
|
|
@ -89,21 +89,8 @@ int uk(
|
|||
// iterate
|
||||
time=starting_time;
|
||||
while(final_time<0 || time<final_time){
|
||||
if(algorithm==ALGORITHM_RK2){
|
||||
ns_step_rk2(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible, keep_en_cst, target_en);
|
||||
} else if (algorithm==ALGORITHM_RK4) {
|
||||
ns_step_rk4(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible, keep_en_cst, target_en);
|
||||
} else if (algorithm==ALGORITHM_RKF45) {
|
||||
ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en, true);
|
||||
} else if (algorithm==ALGORITHM_RKDP54) {
|
||||
// only compute k1 on the first step
|
||||
ns_step_rkdp54(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en, time==starting_time);
|
||||
} else if (algorithm==ALGORITHM_RKBS32) {
|
||||
// only compute k1 on the first step
|
||||
ns_step_rkbs32(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, tmp2, tmp3, &tmp4, tmp5, irreversible, keep_en_cst, target_en, time==starting_time);
|
||||
} else {
|
||||
fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
|
||||
}
|
||||
// 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);
|
||||
|
||||
time+=step;
|
||||
step=next_step;
|
||||
|
@ -205,21 +192,7 @@ int enstrophy(
|
|||
// iterate
|
||||
time=starting_time;
|
||||
while(final_time<0 || time<final_time){
|
||||
if(algorithm==ALGORITHM_RK2){
|
||||
ns_step_rk2(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible, keep_en_cst, target_en);
|
||||
} else if (algorithm==ALGORITHM_RK4) {
|
||||
ns_step_rk4(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible, keep_en_cst, target_en);
|
||||
} else if (algorithm==ALGORITHM_RKF45) {
|
||||
ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en, true);
|
||||
} else if (algorithm==ALGORITHM_RKDP54) {
|
||||
// only compute k1 on the first step
|
||||
ns_step_rkdp54(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en, time==starting_time);
|
||||
} else if (algorithm==ALGORITHM_RKBS32) {
|
||||
// only compute k1 on the first step
|
||||
ns_step_rkbs32(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, tmp2, tmp3, &tmp4, tmp5, irreversible, keep_en_cst, target_en, time==starting_time);
|
||||
} else {
|
||||
fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
|
||||
}
|
||||
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);
|
||||
|
||||
time+=step;
|
||||
|
||||
|
@ -386,21 +359,7 @@ int quiet(
|
|||
// iterate
|
||||
time=starting_time;
|
||||
while(final_time<0 || time<final_time){
|
||||
if(algorithm==ALGORITHM_RK2){
|
||||
ns_step_rk2(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible, keep_en_cst, target_en);
|
||||
} else if (algorithm==ALGORITHM_RK4) {
|
||||
ns_step_rk4(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible, keep_en_cst, target_en);
|
||||
} else if (algorithm==ALGORITHM_RKF45) {
|
||||
ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en, true);
|
||||
} else if (algorithm==ALGORITHM_RKDP54) {
|
||||
// only compute k1 on the first step
|
||||
ns_step_rkdp54(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en, time==starting_time);
|
||||
} else if (algorithm==ALGORITHM_RKBS32) {
|
||||
// only compute k1 on the first step
|
||||
ns_step_rkbs32(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, tmp2, tmp3, &tmp4, tmp5, irreversible, keep_en_cst, target_en, time==starting_time);
|
||||
} else {
|
||||
fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
|
||||
}
|
||||
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);
|
||||
|
||||
time+=step;
|
||||
step=next_step;
|
||||
|
@ -553,6 +512,60 @@ int copy_u(
|
|||
}
|
||||
|
||||
// next time step
|
||||
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_norm,
|
||||
double L,
|
||||
_Complex double* g,
|
||||
double time,
|
||||
double starting_time,
|
||||
fft_vect fft1,
|
||||
fft_vect fft2,
|
||||
fft_vect ifft,
|
||||
// the pointers tmp1 and tmp2 will be exchanged at the end of the routine for first-same-as-last RK algorithms
|
||||
_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
|
||||
){
|
||||
if(algorithm==ALGORITHM_RK2){
|
||||
ns_step_rk2(u, K1, K2, N1, N2, nu, *delta, L, g, fft1, fft2, ifft, *tmp1, *tmp2, irreversible, keep_en_cst, target_en);
|
||||
} else if (algorithm==ALGORITHM_RK4) {
|
||||
ns_step_rk4(u, K1, K2, N1, N2, nu, *delta, L, g, fft1, fft2, ifft, *tmp1, *tmp2, tmp3, irreversible, keep_en_cst, target_en);
|
||||
} else if (algorithm==ALGORITHM_RKF45) {
|
||||
ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, delta, next_delta, L, g, fft1, fft2, ifft, *tmp1, *tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en, true);
|
||||
} else if (algorithm==ALGORITHM_RKDP54) {
|
||||
// only compute k1 on the first step
|
||||
// first-same-as-last with 2-nd argument
|
||||
ns_step_rkdp54(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, delta, next_delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en, time==starting_time);
|
||||
} else if (algorithm==ALGORITHM_RKBS32) {
|
||||
// only compute k1 on the first step
|
||||
// first-same-as-last with 4-th argument
|
||||
ns_step_rkbs32(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, delta, next_delta, L, g, fft1, fft2, ifft, tmp1, tmp3, tmp4, tmp2, tmp5, irreversible, keep_en_cst, target_en, time==starting_time);
|
||||
} else {
|
||||
fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers.edu\n",algorithm);
|
||||
}
|
||||
|
||||
return (0);
|
||||
}
|
||||
|
||||
// RK 4 algorithm
|
||||
int ns_step_rk4(
|
||||
_Complex double* u,
|
||||
|
@ -1342,6 +1355,80 @@ int ns_T_nofft(
|
|||
return 0;
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
// Jacobian of rhs
|
||||
int ns_D_rhs(
|
||||
_Complex double* out,
|
||||
_Complex double* u,
|
||||
int K1,
|
||||
int K2,
|
||||
double nu,
|
||||
double L,
|
||||
bool irreversible
|
||||
){
|
||||
int kx,ky,lx,ly;
|
||||
int i;
|
||||
double alpha;
|
||||
|
||||
if (irreversible) {
|
||||
alpha=nu;
|
||||
} else {
|
||||
// TODO
|
||||
alpha=0;
|
||||
}
|
||||
|
||||
for(i=0; i<MATSIZE*MATSIZE; i++){
|
||||
out[i]=0;
|
||||
}
|
||||
// diagonal term
|
||||
for(kx=0;kx<=K1;kx++){
|
||||
for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
|
||||
out[klookup_sym(kx,ky,K2)*MATSIZE+klookup_sym(kx,ky,K2)]=-4*M_PI*M_PI/L/L*alpha*(kx*kx+ky*ky);
|
||||
}
|
||||
}
|
||||
// convolution term
|
||||
for(kx=0;kx<=K1;kx++){
|
||||
for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
|
||||
for(lx=0;lx<=K1;lx++){
|
||||
for(ly=(lx>0 ? -K2 : 1);ly<=K2;ly++){
|
||||
// no contribution from k=l
|
||||
if (kx!=lx && ky!=ly){
|
||||
// use row-major ordering
|
||||
out[klookup_sym(kx,ky,K2)*MATSIZE+klookup_sym(lx,ly,K2)]=4*M_PI*M_PI/L/L*ns_d_T(u,kx,ky,lx,ly,K2);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return(0);
|
||||
}
|
||||
|
||||
// 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
|
||||
){
|
||||
int qx=kx-lx;
|
||||
int qy=ky-ly;
|
||||
// trivial case
|
||||
if (qx*qx+qy*qy==0){
|
||||
return 0.;
|
||||
}
|
||||
if (qx>0 || (qx==0 && qy>=1)){
|
||||
return (ky*lx-kx*ly)*(qx*qx+qy*qy-lx*lx-ly*ly)/sqrt(lx*lx+ly*ly)/sqrt(qx*qx+qy*qy)*u[klookup_sym(qx,qy,K2)];
|
||||
}else{
|
||||
return (ky*lx-kx*ly)*(qx*qx+qy*qy-lx*lx-ly*ly)/sqrt(lx*lx+ly*ly)/sqrt(qx*qx+qy*qy)*u[klookup_sym(-qx,-qy,K2)];
|
||||
}
|
||||
}
|
||||
*/
|
||||
|
||||
|
||||
// compute alpha
|
||||
double compute_alpha(
|
||||
_Complex double* u,
|
||||
|
|
|
@ -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;
|
||||
|
||||
|
@ -51,6 +49,7 @@ 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_norm, 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
|
||||
|
@ -71,6 +70,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
|
||||
|
|
Loading…
Reference in New Issue