From 527365d62e8c770142c3c1065b5f973433bd60b2 Mon Sep 17 00:00:00 2001
From: Ian Jauslin <ian.jauslin@rutgers.edu>
Date: Wed, 26 Apr 2023 11:13:50 -0400
Subject: Implement RK2

---
 src/navier-stokes.c | 50 ++++++++++++++++++++++++++++++++++++++++++++++----
 1 file changed, 46 insertions(+), 4 deletions(-)

(limited to 'src/navier-stokes.c')

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,
-- 
cgit v1.2.3-54-g00ecf