Implement RK2

2 files changed, 50 insertions, 5 deletions

diff --git a/src/navier-stokes.c b/src/navier-stokes.c
index ce07077..755e829 100644
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -51,7 +51,7 @@ int uk(
 
   // iterate
   for(t=starting_time;nsteps==0 || t<starting_time+nsteps;t++){
-    ns_step(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, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
     
     if(t%print_freq==0){
       fprintf(stderr,"%lu % .8e ",t,t*delta);
@@ -127,7 +127,7 @@ int eea(
 
   // iterate
   for(t=starting_time;nsteps==0 || t<starting_time+nsteps;t++){
-    ns_step(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, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
 
     energy=compute_energy(u, K1, K2);
     alpha=compute_alpha(u, K1, K2, g, L);
@@ -214,7 +214,7 @@ int quiet(
 
   // iterate
   for(t=starting_time;nsteps==0 || t<starting_time+nsteps;t++){
-    ns_step(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, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
   }
 
   // save final entry to savefile
@@ -311,7 +311,8 @@ int copy_u(
 }
 
 // next time step
-int ns_step(
+// RK 4 algorithm
+int ns_step_rk4(
   _Complex double* u,
   int K1,
   int K2,
@@ -388,6 +389,47 @@ int ns_step(
   return(0);
 }
 
+// RK 2 algorithm
+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
+){
+  int kx,ky;
+
+  // k1
+  ns_rhs(tmp1, u, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
+
+  // u+h*k1/2
+  for(kx=0;kx<=K1;kx++){
+    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+      tmp2[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta/2*tmp1[klookup_sym(kx,ky,K2)];
+    }
+  }
+  // k2
+  ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
+  // add to output
+  for(kx=0;kx<=K1;kx++){
+    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+      u[klookup_sym(kx,ky,K2)]+=delta*tmp1[klookup_sym(kx,ky,K2)];
+    }
+  }
+
+  return(0);
+}
+
 // right side of Irreversible/Reversible Navier-Stokes equation
 int ns_rhs(
   _Complex double* out,
diff --git a/src/navier-stokes.h b/src/navier-stokes.h
index d12f908..77b4cdb 100644
--- a/src/navier-stokes.h
+++ b/src/navier-stokes.h
@@ -36,7 +36,10 @@ 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( _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);
+// 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);
+// 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);
 
 // 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);