software/Nstrophy
1
Fork 0
Code Issues Pull Requests Releases Activity

Overhaul of lyapunov exponents

Browse Source
This commit is contained in:
Ian Jauslin 2025-01-31 10:52:40 -05:00
parent fa52397e87
commit 0f07f025b5
5 changed files with 651 additions and 177 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
src/constants.cpp
View File

@ -48,3 +48,6 @@ limitations under the License.
#define FIX_ENSTROPHY 1 #define FIX_ENSTROPHY 1
#define FIX_ENERGY 2 #define FIX_ENERGY 2
#define LYAPUNOV_TRIGGER_TIME 1
#define LYAPUNOV_TRIGGER_SIZE 2

724
src/lyapunov.c
View File

@ -32,110 +32,113 @@ int lyapunov(
int N2, int N2,
double final_time, double final_time,
double lyapunov_reset, double lyapunov_reset,
unsigned int lyapunov_trigger,
double nu, double nu,
double D_epsilon,
double delta, double delta,
double L, double L,
double adaptive_tolerance, double adaptive_tolerance,
double adaptive_factor, double adaptive_factor,
double max_delta, double max_delta,
unsigned int adaptive_norm, unsigned int adaptive_cost,
_Complex double* u0, _Complex double* u0,
_Complex double* g, _Complex double* g,
bool irreversible, bool irreversible,
bool keep_en_cst, bool keep_en_cst,
double target_en, double target_en,
unsigned int algorithm, unsigned int algorithm,
unsigned int algorithm_lyapunov,
double starting_time, double starting_time,
unsigned int nthreads unsigned int nthreads
){ ){
double* lyapunov; double* lyapunov;
double* lyapunov_avg; double* lyapunov_avg;
double* Du_prod; double* flow;
double* Du;
_Complex double* u; _Complex double* u;
_Complex double* prevu; double time;
_Complex double* tmp1; fft_vect fftu1;
_Complex double* tmp2; fft_vect fftu2;
_Complex double* tmp3; fft_vect fftu3;
fft_vect fftu4;
fft_vect fft1;
fft_vect ifft;
double* tmp1;
double* tmp2;
double* tmp3;
_Complex double* tmp4; _Complex double* tmp4;
_Complex double* tmp5; _Complex double* tmp5;
_Complex double* tmp6; _Complex double* tmp6;
_Complex double* tmp7; _Complex double* tmp7;
_Complex double* tmp8; _Complex double* tmp8;
_Complex double* tmp9; _Complex double* tmp9;
double* tmp10; _Complex double* tmp10;
double time; double* tmp11;
fft_vect fft1; int i,j;
fft_vect fft2;
fft_vect ifft;
// period // period (only useful with LYAPUNOV_TRIGGER_TIME, but compute it anyways: it won't take much time...)
// add 0.1 to ensure proper rounding // add 0.1 to ensure proper rounding
uint64_t n=(uint64_t)((starting_time-fmod(starting_time, lyapunov_reset))/lyapunov_reset+0.1); uint64_t n=(uint64_t)((starting_time-fmod(starting_time, lyapunov_reset))/lyapunov_reset+0.1);
lyapunov_init_tmps(&lyapunov, &lyapunov_avg, &Du_prod, &Du, &u, &prevu, &tmp1, &tmp2, &tmp3, &tmp4, &tmp5, &tmp6, &tmp7, &tmp8, &tmp9, &tmp10, &fft1, &fft2, &ifft, K1, K2, N1, N2, nthreads, algorithm); lyapunov_init_tmps(&lyapunov, &lyapunov_avg, &flow, &u, &tmp1, &tmp2, &tmp3, &tmp4, &tmp5, &tmp6, &tmp7, &tmp8, &tmp9, &tmp10, &tmp11, &fftu1, &fftu2, &fftu3, &fftu4, &fft1, &ifft, K1, K2, N1, N2, nthreads, algorithm, algorithm_lyapunov);
// copy initial condition // copy initial condition
copy_u(u, u0, K1, K2); copy_u(u, u0, K1, K2);
// initialize flow with identity matrix
// initialize Du_prod for (i=0;i<MATSIZE;i++){
int i,j; for (j=0;j<MATSIZE;j++){
for(i=0;i<MATSIZE;i++){ if(i!=j){
for(j=0;j<MATSIZE;j++){ flow[i*MATSIZE+j]=0.;
Du_prod[i*MATSIZE+j]=0.; } else {
flow[i*MATSIZE+j]=1.;
} }
} }
for(i=0;i<MATSIZE;i++){
Du_prod[i*MATSIZE+i]=1.;
} }
// store step (useful for adaptive step methods // store step (useful for adaptive step methods)
double prev_step=delta;
double step=delta; double step=delta;
double next_step=step; double next_step=step;
// init average
for (i=0; i<MATSIZE; i++){
lyapunov_avg[i]=0;
}
// save times at which Lyapunov exponents are computed // save times at which Lyapunov exponents are computed
double prevtime=starting_time; double prevtime=starting_time;
// iterate // iterate
time=starting_time; time=starting_time;
while(final_time<0 || time<final_time){ while(final_time<0 || time<final_time){
// copy before step // compute u first
copy_u(prevu, u, K1, K2); // if using an adaptive step, the step for the tangent flow is set by this computation of u
prev_step=step; // 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);
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 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);
// compute Jacobian
// do not touch tmp1, it might be used in the next step
ns_D_step(Du, prevu, u, K1, K2, N1, N2, nu, D_epsilon, prev_step, algorithm, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, L, g, time, fft1, fft2, ifft, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9, irreversible, keep_en_cst, target_en);
// multiply Jacobian
// switch to column major order, so transpose Du
cblas_dgemm(CblasColMajor, CblasNoTrans, CblasTrans, MATSIZE, MATSIZE, MATSIZE, 1., Du_prod, MATSIZE, Du, MATSIZE, 0., tmp10, MATSIZE);
// copy back to Du_prod
double* move=tmp10;
tmp10=Du_prod;
Du_prod=move;
// increment time // increment time
time+=step; time+=step;
// reset Jacobian double norm=0;
if(time>(n+1)*lyapunov_reset){ if(lyapunov_trigger==LYAPUNOV_TRIGGER_SIZE){
// size of flow (for reset)
for(i=0;i<MATSIZE;i++){
for(j=0;j<MATSIZE;j++){
norm+=fabs(flow[i*MATSIZE+j]);
if(norm>lyapunov_reset){
break;
}
}
if(norm>lyapunov_reset){
break;
}
}
}
// QR decomposition
if((lyapunov_trigger==LYAPUNOV_TRIGGER_TIME && time>(n+1)*lyapunov_reset) || (lyapunov_trigger==LYAPUNOV_TRIGGER_SIZE && norm>lyapunov_reset)){
n++; n++;
// QR decomposition // QR decomposition
// do not touch tmp1, it might be used in the next step // do not touch tmp1, it might be used in the next step
LAPACKE_dgerqf(LAPACK_COL_MAJOR, MATSIZE, MATSIZE, Du_prod, MATSIZE, tmp10); LAPACKE_dgeqrf(LAPACK_COL_MAJOR, MATSIZE, MATSIZE, flow, MATSIZE, tmp11);
// extract eigenvalues (diagonal elements of Du_prod) // extract diagonal elements of R (represented as diagonal elements of flow
for(i=0; i<MATSIZE; i++){ for(i=0; i<MATSIZE; i++){
lyapunov[i]=log(fabs(Du_prod[i*MATSIZE+i]))/(time-prevtime); lyapunov[i]=log(fabs(flow[i*MATSIZE+i]))/(time-prevtime);
} }
// sort lyapunov exponents // sort lyapunov exponents
@ -155,23 +158,24 @@ int lyapunov(
} }
printf("\n"); printf("\n");
fprintf(stderr,"% .15e",time); fprintf(stderr,"% .15e",time);
// print largest and smallest lyapunov exponent // print largest and smallest lyapunov exponent to stderr
if(MATSIZE>0){ if(MATSIZE>0){
fprintf(stderr," % .15e % .15e\n", lyapunov[0], lyapunov[MATSIZE-1]); fprintf(stderr," % .15e % .15e\n", lyapunov[0], lyapunov[MATSIZE-1]);
} }
// set Du_prod to Q: // set flow to Q:
LAPACKE_dorgrq(LAPACK_COL_MAJOR, MATSIZE, MATSIZE, MATSIZE, Du_prod, MATSIZE, tmp10); LAPACKE_dorgqr(LAPACK_COL_MAJOR, MATSIZE, MATSIZE, MATSIZE, flow, MATSIZE, tmp11);
// reset prevtime // reset prevtime
prevtime=time; 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); 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); return(0);
} }
// comparison function for qsort // comparison function for qsort
int compare_double(const void* x, const void* y) { int compare_double(const void* x, const void* y) {
if (*(const double*)x<*(const double*)y) { if (*(const double*)x<*(const double*)y) {
@ -183,87 +187,64 @@ int compare_double(const void* x, const void* y) {
} }
} }
// Jacobian of u_n -> u_{n+1} // compute tangent flow
int ns_D_step( int lyapunov_D_step(
double* out, double* flow,
_Complex double* un, _Complex double* u,
_Complex double* un_next,
int K1, int K1,
int K2, int K2,
int N1, int N1,
int N2, int N2,
double nu, double nu,
double epsilon,
double delta, double delta,
unsigned int algorithm, unsigned int algorithm,
double adaptive_tolerance,
double adaptive_factor,
double max_delta,
unsigned int adaptive_norm,
double L, double L,
_Complex double* g, _Complex double* g,
double time, fft_vect fftu1,
fft_vect fftu2,
fft_vect fftu3,
fft_vect fftu4,
fft_vect fft1, fft_vect fft1,
fft_vect fft2,
fft_vect ifft, fft_vect ifft,
_Complex double* tmp1, double* tmp1,
_Complex double* tmp2, double* tmp2,
_Complex double* tmp3, double* tmp3,
_Complex double* tmp4, _Complex double* tmp4,
_Complex double* tmp5, bool irreversible
_Complex double* tmp6,
_Complex double* tmp7,
_Complex double* tmp8,
bool irreversible,
bool keep_en_cst,
double target_en
){ ){
int lx,ly; int lx,ly;
int i; double alpha;
double step, next_step;
// 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=-K1;lx<=K1;lx++){
for(ly=-K2;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(lx=0;lx<=K1;lx++){
for(ly=(lx>0 ? -K2 : 1);ly<=K2;ly++){ for(ly=(lx>0 ? -K2 : 1);ly<=K2;ly++){
// derivative in the real direction
// finite difference vector
for (i=0;i<USIZE;i++){
if(i==klookup_sym(lx,ly,K2)){
tmp1[i]=un[i]+epsilon;
}else{
tmp1[i]=un[i];
}
}
// compute step
// reinitialize step
step=delta;
next_step=delta;
ns_step(algorithm, tmp1, K1, K2, N1, N2, nu, &step, &next_step, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, L, g, time, time, fft1, fft2, ifft, &tmp2, &tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, irreversible, keep_en_cst, target_en);
// derivative
for (i=0;i<USIZE;i++){
// use row major order
out[2*klookup_sym(lx,ly,K2)*MATSIZE+2*i]=(__real__ (tmp1[i]-un_next[i]))/epsilon;
out[2*klookup_sym(lx,ly,K2)*MATSIZE+2*i+1]=(__imag__ (tmp1[i]-un_next[i]))/epsilon;
}
// derivative in the imaginary direction // do not use adaptive step algorithms for the tangent flow: it must be locked to the same times as u
// finite difference vector if(algorithm==ALGORITHM_RK2){
for (i=0;i<USIZE;i++){ lyapunov_D_step_rk2(flow+2*klookup_sym(lx,ly,K2)*MATSIZE, u, K1, K2, N1, N2, alpha, delta, L, g, tmp4, fftu1, fftu2, fftu3, fftu4, fft1, ifft, tmp1, tmp2, irreversible);
if(i==klookup_sym(lx,ly,K2)){ } else if (algorithm==ALGORITHM_RK4) {
tmp1[i]=un[i]+epsilon*I; lyapunov_D_step_rk4(flow+2*klookup_sym(lx,ly,K2)*MATSIZE, u, K1, K2, N1, N2, alpha, delta, L, g, tmp4, fftu1, fftu2, fftu3, fftu4, fft1, ifft, tmp1, tmp2, tmp3, irreversible);
}else{ } else {
tmp1[i]=un[i]; fprintf(stderr,"bug: unknown algorithm for tangent flow: %u, contact ian.jauslin@rutgers.edu\n",algorithm);
}
}
// compute step
// reinitialize step
step=delta;
next_step=delta;
ns_step(algorithm, tmp1, K1, K2, N1, N2, nu, &step, &next_step, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, L, g, time, time, fft1, fft2, ifft, &tmp2, &tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, irreversible, keep_en_cst, target_en);
// compute derivative
for (i=0;i<USIZE;i++){
// use row major order
out[(2*klookup_sym(lx,ly,K2)+1)*MATSIZE+2*i]=(__real__ (tmp1[i]-un_next[i]))/epsilon;
out[(2*klookup_sym(lx,ly,K2)+1)*MATSIZE+2*i+1]=(__imag__ (tmp1[i]-un_next[i]))/epsilon;
} }
} }
} }
@ -271,80 +252,501 @@ int ns_D_step(
return(0); 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( int lyapunov_init_tmps(
double ** lyapunov, double ** lyapunov,
double ** lyapunov_avg, double ** lyapunov_avg,
double ** Du_prod, double ** flow,
double ** Du,
_Complex double ** u, _Complex double ** u,
_Complex double ** prevu, double ** tmp1,
_Complex double ** tmp1, double ** tmp2,
_Complex double ** tmp2, double ** tmp3,
_Complex double ** tmp3,
_Complex double ** tmp4, _Complex double ** tmp4,
_Complex double ** tmp5, _Complex double ** tmp5,
_Complex double ** tmp6, _Complex double ** tmp6,
_Complex double ** tmp7, _Complex double ** tmp7,
_Complex double ** tmp8, _Complex double ** tmp8,
_Complex double ** tmp9, _Complex double ** tmp9,
double ** tmp10, _Complex double ** tmp10,
double ** tmp11,
fft_vect* fftu1,
fft_vect* fftu2,
fft_vect* fftu3,
fft_vect* fftu4,
fft_vect* fft1, fft_vect* fft1,
fft_vect* fft2,
fft_vect* ifft, fft_vect* ifft,
int K1, int K1,
int K2, int K2,
int N1, int N1,
int N2, int N2,
unsigned int nthreads, unsigned int nthreads,
unsigned int algorithm unsigned int algorithm,
unsigned int algorithm_lyapunov
){ ){
ns_init_tmps(u, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9, fft1, fft2, ifft, K1, K2, N1, N2, nthreads, algorithm); // velocity field
*u=calloc(USIZE,sizeof(_Complex double));
// 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=calloc(MATSIZE,sizeof(double));
*lyapunov_avg=calloc(MATSIZE,sizeof(double)); *lyapunov_avg=calloc(MATSIZE,sizeof(double));
*Du_prod=calloc(MATSIZE*MATSIZE,sizeof(double)); *flow=calloc(MATSIZE*MATSIZE,sizeof(double));
*Du=calloc(MATSIZE*MATSIZE,sizeof(double));
*prevu=calloc(USIZE,sizeof(_Complex double));
*tmp1=calloc(USIZE,sizeof(_Complex double));
*tmp2=calloc(USIZE,sizeof(_Complex double));
*tmp10=calloc(MATSIZE*MATSIZE,sizeof(double));
return(0); return(0);
} }
// release vectors // release vectors
int lyapunov_free_tmps( int lyapunov_free_tmps(
double* lyapunov, double * lyapunov,
double* lyapunov_avg, double * lyapunov_avg,
double* Du_prod, double * flow,
double* Du, _Complex double * u,
_Complex double* u, double * tmp1,
_Complex double* prevu, double * tmp2,
_Complex double* tmp1, double * tmp3,
_Complex double* tmp2, _Complex double * tmp4,
_Complex double* tmp3, _Complex double * tmp5,
_Complex double* tmp4, _Complex double * tmp6,
_Complex double* tmp5, _Complex double * tmp7,
_Complex double* tmp6, _Complex double * tmp8,
_Complex double* tmp7, _Complex double * tmp9,
_Complex double* tmp8, _Complex double * tmp10,
_Complex double* tmp9, double * tmp11,
double* tmp10, fft_vect fftu1,
fft_vect fftu2,
fft_vect fftu3,
fft_vect fftu4,
fft_vect fft1, fft_vect fft1,
fft_vect fft2,
fft_vect ifft, fft_vect ifft,
unsigned int algorithm unsigned int algorithm,
unsigned int algorithm_lyapunov
){ ){
free(tmp10); free(flow);
free(tmp2);
free(tmp1);
free(prevu);
free(Du);
free(Du_prod);
free(lyapunov_avg);
free(lyapunov); free(lyapunov);
ns_free_tmps(u, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9, fft1, fft2, ifft, algorithm); free(lyapunov_avg);
fftw_destroy_plan(fftu2.fft_plan);
fftw_destroy_plan(fftu3.fft_plan);
fftw_destroy_plan(fftu4.fft_plan);
fftw_free(fftu2.fft);
fftw_free(fftu3.fft);
fftw_free(fftu4.fft);
ns_free_tmps(u, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9, tmp10, fft1, fftu1, ifft, algorithm);
free(tmp11);
if(algorithm_lyapunov==ALGORITHM_RK2){
free(tmp1);
free(tmp2);
} else if (algorithm_lyapunov==ALGORITHM_RK4){
free(tmp1);
free(tmp2);
free(tmp3);
} else {
fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm_lyapunov);
};
return(0); return(0);
} }

35
src/lyapunov.h
View File

@ -20,18 +20,41 @@ limitations under the License.
#include "navier-stokes.h" #include "navier-stokes.h"
// compute Lyapunov exponents // compute Lyapunov exponents
int lyapunov( int K1, int K2, int N1, int N2, double final_time, double lyapunov_reset, double nu, double D_epsilon, double delta, double L, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_norm, _Complex double* u0, _Complex double* g, bool irreversible, bool keep_en_cst, double target_en, unsigned int algorithm, double starting_time, unsigned int nthreads); int lyapunov( int K1, int K2, int N1, int N2, double final_time, double lyapunov_reset, unsigned int lyapunov_trigger, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_norm, _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);
// comparison function for qsort // comparison function for qsort
int compare_double(const void* x, const void* y); int compare_double(const void* x, const void* y);
// Jacobian of step // compute tangent flow
int ns_D_step( double* out, _Complex double* un, _Complex double* un_next, int K1, int K2, int N1, int N2, double nu, double epsilon, double delta, unsigned int algorithm, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_norm, double L, _Complex double* g, double time, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double* tmp1, _Complex double* tmp2, _Complex double* tmp3, _Complex double* tmp4, _Complex double* tmp5, _Complex double* tmp6, _Complex double* tmp7, _Complex double* tmp8, bool irreversible, bool keep_en_cst, double target_en); int lyapunov_D_step( double* flow, _Complex double* u, int K1, int K2, int N1, int N2, double nu, double delta, unsigned int algorithm, double L, _Complex double* g, fft_vect fftu1, fft_vect fftu2, fft_vect fftu3, fft_vect fftu4, fft_vect fft1, fft_vect ifft, double* tmp1, double* tmp2, double* tmp3, _Complex double* tmp4, bool irreversible);
// RK 4 algorithm
int lyapunov_D_step_rk4( double* flow, _Complex double* u, int K1, int K2, int N1, int N2, double nu, double delta, double L, _Complex double* g, _Complex double* T, fft_vect fftu1, fft_vect fftu2, fft_vect fftu3, fft_vect fftu4, fft_vect fft1, fft_vect ifft, double* tmp1, double* tmp2, double* tmp3, bool irreversible);
// RK 2 algorithm
int lyapunov_D_step_rk2( double* flow, _Complex double* u, int K1, int K2, int N1, int N2, double nu, double delta, double L, _Complex double* g, _Complex double* T, fft_vect fftu1, fft_vect fftu2, fft_vect fftu3, fft_vect fftu4, fft_vect fft1, fft_vect ifft, double* tmp1, double* tmp2, bool irreversible);
// Right side of equation for tangent flow
int lyapunov_D_rhs( double* out, double* flow, _Complex double* u, int K1, int K2, int N1, int N2, double alpha, double L, _Complex double* g, _Complex double* T, fft_vect fftu1, fft_vect fftu2, fft_vect fftu3, fft_vect fftu4, fft_vect fft1, fft_vect ifft, bool irreversible);
// Compute T from the already computed fourier transforms of u
int lyapunov_T( int N1, int N2, fft_vect fftu1, fft_vect fftu2, fft_vect fftu3, fft_vect fftu4, fft_vect ifft);
// compute derivative of T
int lyapunov_D_T( double* flow, int K1, int K2, int N1, int N2, fft_vect fftu1, fft_vect fftu2, fft_vect fftu3, fft_vect fftu4, fft_vect fft1, fft_vect ifft);
double lyapunov_D_alpha( double* flow, _Complex double* u, int K1, int K2, int N1, int N2, double alpha, double L, _Complex double* g, _Complex double* T, _Complex double* DT);
// get delta_{kx,ky} from a vector delta in which only the values for kx>=0 are stored, as separate real part and imaginary part
_Complex double lyapunov_delta_getval_sym( double* delta, int kx, int ky, int K2);
// compute ffts of u for future use in DT
int lyapunov_fft_u_T( _Complex double* u, int K1, int K2, int N1, int N2, fft_vect fftu1, fft_vect fftu2, fft_vect fftu3, fft_vect fftu4);
int lyapunov_init_tmps( double ** lyapunov, double ** lyapunov_avg, double ** flow, _Complex double ** u, double ** tmp1, double ** tmp2, double ** tmp3, _Complex double ** tmp4, _Complex double ** tmp5, _Complex double ** tmp6, _Complex double ** tmp7, _Complex double ** tmp8, _Complex double ** tmp9, _Complex double ** tmp10, double ** tmp11, fft_vect* fftu1, fft_vect* fftu2, fft_vect* fftu3, fft_vect* fftu4, fft_vect* fft1, fft_vect* ifft, int K1, int K2, int N1, int N2, unsigned int nthreads, unsigned int algorithm, unsigned int algorithm_lyapunov);
// init vectors
int lyapunov_init_tmps(double** lyapunov, double** lyapunov_avg, double ** Du_prod, double ** Du, _Complex double ** u, _Complex double ** prevu, _Complex double ** tmp1, _Complex double ** tmp2, _Complex double ** tmp3, _Complex double ** tmp4, _Complex double ** tmp5, _Complex double ** tmp6, _Complex double ** tmp7, _Complex double ** tmp8, _Complex double ** tmp9, double** tmp10, fft_vect* fft1, fft_vect* fft2, fft_vect* ifft, int K1, int K2, int N1, int N2, unsigned int nthreads, unsigned int algorithm);
// release vectors // release vectors
int lyapunov_free_tmps(double* lyapunov, double* lyapunov_avg, double* Du_prod, double* Du, _Complex double* u, _Complex double* prevu, _Complex double* tmp1, _Complex double* tmp2, _Complex double* tmp3, _Complex double* tmp4, _Complex double* tmp5, _Complex double* tmp6, _Complex double* tmp7, _Complex double* tmp8, _Complex double* tmp9, double* tmp10, fft_vect fft1, fft_vect fft2, fft_vect ifft, unsigned int algorithm); int lyapunov_free_tmps( double * lyapunov, double * lyapunov_avg, double * flow, _Complex double * u, double * tmp1, double * tmp2, double * tmp3, _Complex double * tmp4, _Complex double * tmp5, _Complex double * tmp6, _Complex double * tmp7, _Complex double * tmp8, _Complex double * tmp9, _Complex double * tmp10, double * tmp11, fft_vect fftu1, fft_vect fftu2, fft_vect fftu3, fft_vect fftu4, fft_vect fft1, fft_vect ifft, unsigned int algorithm, unsigned int algorithm_lyapunov);
#endif #endif

61
src/main.c
View File

@ -56,10 +56,12 @@ typedef struct nstrophy_parameters {
unsigned int driving; unsigned int driving;
unsigned int init; unsigned int init;
unsigned int algorithm; unsigned int algorithm;
unsigned int algorithm_lyapunov;
bool keep_en_cst; bool keep_en_cst;
FILE* initfile; FILE* initfile;
FILE* drivingfile; FILE* drivingfile;
double lyapunov_reset; double lyapunov_reset;
unsigned int lyapunov_trigger;
double D_epsilon; double D_epsilon;
bool print_alpha; bool print_alpha;
} nstrophy_parameters; } nstrophy_parameters;
@ -296,7 +298,7 @@ int main (
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); 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){ else if(command==COMMAND_LYAPUNOV){
lyapunov(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.lyapunov_reset, parameters.nu, parameters.D_epsilon, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_cost, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_enstrophy, parameters.algorithm, parameters.starting_time, nthreads); 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);
} }
else if(command==0){ else if(command==0){
fprintf(stderr, "error: no command specified\n"); fprintf(stderr, "error: no command specified\n");
@ -576,6 +578,9 @@ int set_default_params(
parameters->initfile=NULL; parameters->initfile=NULL;
parameters->algorithm=ALGORITHM_RK4; parameters->algorithm=ALGORITHM_RK4;
parameters->keep_en_cst=false; parameters->keep_en_cst=false;
parameters->algorithm_lyapunov=ALGORITHM_RK4;
// default lyapunov_reset will be print_time (set later) for now set target to 0 to indicate it hasn't been set explicitly
parameters->lyapunov_trigger=0;
parameters->print_alpha=false; parameters->print_alpha=false;
@ -648,6 +653,12 @@ int read_params(
parameters->N2=smallest_pow2(3*(parameters->K2)); 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); return(0);
} }
@ -841,13 +852,6 @@ int set_parameter(
return(-1); 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){ else if (strcmp(lhs,"D_epsilon")==0){
ret=sscanf(rhs,"%lf",&(parameters->D_epsilon)); ret=sscanf(rhs,"%lf",&(parameters->D_epsilon));
if(ret!=1){ if(ret!=1){
@ -940,6 +944,47 @@ int set_parameter(
} }
parameters->print_alpha=(tmp==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{ else{
fprintf(stderr, "error: unrecognized parameter '%s'\n",lhs); fprintf(stderr, "error: unrecognized parameter '%s'\n",lhs);
return(-1); return(-1);

3
src/navier-stokes.c
View File

@ -82,7 +82,7 @@ int uk(
// add 0.1 to ensure proper rounding // add 0.1 to ensure proper rounding
uint64_t n=(uint64_t)((starting_time-fmod(starting_time, print_freq))/print_freq+0.1); uint64_t n=(uint64_t)((starting_time-fmod(starting_time, print_freq))/print_freq+0.1);
// store step (useful for adaptive step methods // store step (useful for adaptive step methods)
double step=delta; double step=delta;
double next_step=step; double next_step=step;
@ -524,6 +524,7 @@ int ns_step(
return (0); return (0);
} }
// TODO: do not need to use klookup in any of the rk routines
// RK 4 algorithm // RK 4 algorithm
int ns_step_rk4( int ns_step_rk4(
_Complex double* u, _Complex double* u,