From dacd5c36511c95fb22658f1634e49ceb1e5c6fd2 Mon Sep 17 00:00:00 2001
From: Ian Jauslin <ian.jauslin@rutgers.edu>
Date: Tue, 16 May 2023 17:03:58 -0400
Subject: [PATCH] Adaptive timestep: do nto change delta before all the
 printing

---
 src/navier-stokes.c | 85 ++++++++++++++++++++++++++++-----------------
 src/navier-stokes.h |  2 +-
 2 files changed, 54 insertions(+), 33 deletions(-)

diff --git a/src/navier-stokes.c b/src/navier-stokes.c
index 8cfaa7a..af780a2 100644
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -77,20 +77,25 @@ int uk(
   // add 0.1 to ensure proper rounding
   uint64_t n=(uint64_t)((starting_time-fmod(starting_time, print_freq))/print_freq+0.1);
 
+  // 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){
     if(algorithm==ALGORITHM_RK2){
-      ns_step_rk2(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible);
+      ns_step_rk2(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible);
     } else if (algorithm==ALGORITHM_RK4) {
-      ns_step_rk4(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
+      ns_step_rk4(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
     } else if (algorithm==ALGORITHM_RKF45) {
-      delta=ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible);
+      ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible);
     } else {
       fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
     }
 
-    time+=delta;
+    time+=step;
+    step=next_step;
     
     if(time>(n+1)*print_freq){
       n++;
@@ -165,6 +170,9 @@ int enstrophy(
   // copy initial condition
   copy_u(u, u0, K1, K2);
 
+  // store step (useful for adaptive step methods
+  double step=delta;
+  double next_step=step;
 
   // init running average
   avg_a=0;
@@ -180,24 +188,24 @@ int enstrophy(
   time=starting_time;
   while(final_time<0 || time<final_time){
     if(algorithm==ALGORITHM_RK2){
-      ns_step_rk2(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible);
+      ns_step_rk2(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible);
     } else if (algorithm==ALGORITHM_RK4) {
-      ns_step_rk4(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
+      ns_step_rk4(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
     } else if (algorithm==ALGORITHM_RKF45) {
-      delta=ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible);
+      ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible);
     } else {
       fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
     }
 
-    time+=delta;
+    time+=step;
 
     alpha=compute_alpha(u, K1, K2, g, L);
     enstrophy=compute_enstrophy(u, K1, K2, L);
 
     // add to running averages (estimate the total duration of interval as print_freq, will be adjusted later)
-    avg_a+=alpha*(delta/print_freq);
-    avg_en+=enstrophy*(delta/print_freq);
-    avg_en_x_a+=enstrophy*alpha*(delta/print_freq);
+    avg_a+=alpha*(step/print_freq);
+    avg_en+=enstrophy*(step/print_freq);
+    avg_en_x_a+=enstrophy*alpha*(step/print_freq);
 
     if(time>(n+1)*print_freq){
       n++;
@@ -209,8 +217,8 @@ int enstrophy(
 
       // print to stderr so user can follow along
       if(algorithm==ALGORITHM_RKF45){
-	fprintf(stderr,"% .8e % .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",time, avg_a, avg_en, avg_en_x_a, alpha, enstrophy, alpha*enstrophy, delta);
-	printf("% .15e % .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",time, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy, delta);
+	fprintf(stderr,"% .8e % .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",time, avg_a, avg_en, avg_en_x_a, alpha, enstrophy, alpha*enstrophy, step);
+	printf("% .15e % .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",time, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy, step);
       } else {
 	fprintf(stderr,"% .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",time, avg_a, avg_en, avg_en_x_a, alpha, enstrophy, alpha*enstrophy);
 	printf("% .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",time, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy);
@@ -224,6 +232,10 @@ int enstrophy(
     avg_en_x_a=0;
     prevtime=time;
 
+    // get ready for next step
+    step=next_step;
+
+
     // catch abort signal
     if (g_abort){
       // print u to stderr if no savefile
@@ -254,7 +266,7 @@ int enstrophy(
       free(params);
       // write delta if adaptive, and not writing binary
       if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD && (savefile==stderr || savefile==stdout)){
-	fprintf(savefile,";delta=%.15e", delta);
+	fprintf(savefile,";delta=%.15e", step);
       }
       // write starting_time if not writing binary
       if(savefile==stderr || savefile==stdout){
@@ -273,7 +285,7 @@ int enstrophy(
       fwrite(&time, sizeof(double), 1, savefile);
       // extra binary data for adaptive algorithm
       if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){
-	fwrite(&delta, sizeof(double), 1, savefile);
+	fwrite(&step, sizeof(double), 1, savefile);
       }
     }
   }
@@ -319,20 +331,26 @@ int quiet(
   // copy initial condition
   copy_u(u, u0, K1, K2);
 
+  // store delta (useful for adaptive step algorithms)
+  double step=delta;
+  double next_step=step;
+
+
   // iterate
   time=starting_time;
   while(final_time<0 || time<final_time){
     if(algorithm==ALGORITHM_RK2){
-      ns_step_rk2(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible);
+      ns_step_rk2(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible);
     } else if (algorithm==ALGORITHM_RK4) {
-      ns_step_rk4(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
+      ns_step_rk4(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
     } else if (algorithm==ALGORITHM_RKF45) {
-      delta=ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible);
+      ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible);
     } else {
       fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
     }
 
-    time+=delta;
+    time+=step;
+    step=next_step;
   }
 
   // save final entry to savefile
@@ -590,7 +608,7 @@ int ns_step_rk2(
 
 // next time step
 // adaptive RK algorithm (Runge-Kutta-Fehlberg)
-double ns_step_rkf45(
+int ns_step_rkf45(
   _Complex double* u,
   double tolerance,
   double factor,
@@ -599,7 +617,8 @@ double ns_step_rkf45(
   int N1,
   int N2,
   double nu,
-  double delta,
+  double* delta,
+  double* next_delta,
   double L,
   _Complex double* g,
   fft_vect fft1,
@@ -623,7 +642,7 @@ double ns_step_rkf45(
   // k2 : u(t+1/4*delta)
   for(kx=0;kx<=K1;kx++){
     for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta/4*k1[klookup_sym(kx,ky,K2)];
+      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)/4*k1[klookup_sym(kx,ky,K2)];
     }
   }
   ns_rhs(k2, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
@@ -631,7 +650,7 @@ double ns_step_rkf45(
   // k3 : u(t+3/8*delta)
   for(kx=0;kx<=K1;kx++){
     for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta*(3./32*k1[klookup_sym(kx,ky,K2)]+9./32*k2[klookup_sym(kx,ky,K2)]);
+      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(3./32*k1[klookup_sym(kx,ky,K2)]+9./32*k2[klookup_sym(kx,ky,K2)]);
     }
   }
   ns_rhs(k3, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
@@ -639,7 +658,7 @@ double ns_step_rkf45(
   // k4 : u(t+12./13*delta)
   for(kx=0;kx<=K1;kx++){
     for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta*(1932./2197*k1[klookup_sym(kx,ky,K2)]-7200./2197*k2[klookup_sym(kx,ky,K2)]+7296./2197*k3[klookup_sym(kx,ky,K2)]);
+      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(1932./2197*k1[klookup_sym(kx,ky,K2)]-7200./2197*k2[klookup_sym(kx,ky,K2)]+7296./2197*k3[klookup_sym(kx,ky,K2)]);
     }
   }
   ns_rhs(k4, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
@@ -647,7 +666,7 @@ double ns_step_rkf45(
   // k5 : u(t+1*delta)
   for(kx=0;kx<=K1;kx++){
     for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta*(439./216*k1[klookup_sym(kx,ky,K2)]-8*k2[klookup_sym(kx,ky,K2)]+3680./513*k3[klookup_sym(kx,ky,K2)]-845./4104*k4[klookup_sym(kx,ky,K2)]);
+      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(439./216*k1[klookup_sym(kx,ky,K2)]-8*k2[klookup_sym(kx,ky,K2)]+3680./513*k3[klookup_sym(kx,ky,K2)]-845./4104*k4[klookup_sym(kx,ky,K2)]);
     }
   }
   ns_rhs(k5, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
@@ -655,7 +674,7 @@ double ns_step_rkf45(
   // k6 : u(t+1./2*delta)
   for(kx=0;kx<=K1;kx++){
     for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta*(-8./27*k1[klookup_sym(kx,ky,K2)]+2*k2[klookup_sym(kx,ky,K2)]-3544./2565*k3[klookup_sym(kx,ky,K2)]+1859./4104*k4[klookup_sym(kx,ky,K2)]-11./40*k5[klookup_sym(kx,ky,K2)]);
+      tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(-8./27*k1[klookup_sym(kx,ky,K2)]+2*k2[klookup_sym(kx,ky,K2)]-3544./2565*k3[klookup_sym(kx,ky,K2)]+1859./4104*k4[klookup_sym(kx,ky,K2)]-11./40*k5[klookup_sym(kx,ky,K2)]);
     }
   }
   ns_rhs(k6, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
@@ -666,15 +685,13 @@ double ns_step_rkf45(
   for(kx=0;kx<=K1;kx++){
     for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
       // difference between 5th order and 4th order
-      err+=cabs(delta*(1./360*k1[klookup_sym(kx,ky,K2)]-128./4275*k3[klookup_sym(kx,ky,K2)]-2197./75240*k4[klookup_sym(kx,ky,K2)]+1./50*k5[klookup_sym(kx,ky,K2)]+2./55*k6[klookup_sym(kx,ky,K2)]));
+      err+=cabs((*delta)*(1./360*k1[klookup_sym(kx,ky,K2)]-128./4275*k3[klookup_sym(kx,ky,K2)]-2197./75240*k4[klookup_sym(kx,ky,K2)]+1./50*k5[klookup_sym(kx,ky,K2)]+2./55*k6[klookup_sym(kx,ky,K2)]));
       // next step
-      tmp[klookup_sym(kx,ky,K2)]=delta*(16./135*k1[klookup_sym(kx,ky,K2)]+6656./12825*k3[klookup_sym(kx,ky,K2)]+28561./56430*k4[klookup_sym(kx,ky,K2)]-9./50*k5[klookup_sym(kx,ky,K2)]+2./55*k6[klookup_sym(kx,ky,K2)]);
+      tmp[klookup_sym(kx,ky,K2)]=(*delta)*(16./135*k1[klookup_sym(kx,ky,K2)]+6656./12825*k3[klookup_sym(kx,ky,K2)]+28561./56430*k4[klookup_sym(kx,ky,K2)]-9./50*k5[klookup_sym(kx,ky,K2)]+2./55*k6[klookup_sym(kx,ky,K2)]);
       relative+=cabs(tmp[klookup_sym(kx,ky,K2)]);
     }
   }
 
-  // new delta
-  delta=factor*delta*pow(relative*tolerance/err,0.2);
   // compare relative error with tolerance
   if(err<relative*tolerance){
     // add to output
@@ -683,12 +700,16 @@ double ns_step_rkf45(
 	u[klookup_sym(kx,ky,K2)]+=tmp[klookup_sym(kx,ky,K2)];
       }
     }
-    return delta;
+    // next delta to use in future steps
+    *next_delta=factor*(*delta)*pow(relative*tolerance/err,0.2);
   }
   // error too big: repeat with smaller step
   else{
-    return(ns_step_rkf45(u,tolerance,factor,K1,K2,N1,N2,nu,delta,L,g,fft1,fft2,ifft,k1,k2,k3,k4,k5,k6,tmp,irreversible));
+    *delta=factor*(*delta)*pow(relative*tolerance/err,0.2);
+    ns_step_rkf45(u,tolerance,factor,K1,K2,N1,N2,nu,delta,next_delta,L,g,fft1,fft2,ifft,k1,k2,k3,k4,k5,k6,tmp,irreversible);
   }
+
+  return 0;
 }
 
 // right side of Irreversible/Reversible Navier-Stokes equation
diff --git a/src/navier-stokes.h b/src/navier-stokes.h
index 074e744..eff95fb 100644
--- a/src/navier-stokes.h
+++ b/src/navier-stokes.h
@@ -56,7 +56,7 @@ int ns_step_rk4( _Complex double* u, int K1, int K2, int N1, int N2, double nu,
 // 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);
 // adaptive RK algorithm (Runge-Kutta-Fehlberg)
-double ns_step_rkf45( _Complex double* u, double tolerance, double factor, 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* k1, _Complex double* k2, _Complex double* k3, _Complex double* k4, _Complex double* k5, _Complex double* k6, _Complex double* tmp, bool irreversible);
+int ns_step_rkf45( _Complex double* u, double tolerance, double factor, 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);
 
 // 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);