Reversible equation

1 files changed, 60 insertions, 17 deletions

diff --git a/src/navier-stokes.c b/src/navier-stokes.c
index 1980207..07eb343 100644
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -15,6 +15,7 @@ int uk(
   double L,
   _Complex double* u0,
   _Complex double* g,
+  bool irreversible,
   unsigned int print_freq,
   unsigned int starting_time,
   unsigned int nthreads,
@@ -48,7 +49,7 @@ int uk(
 
   // iterate
   for(t=starting_time;t<starting_time+nsteps;t++){
-    ins_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
+    ns_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
     
     if(t%print_freq==0){
       fprintf(stderr,"%d % .8e ",t,t*delta);
@@ -87,6 +88,7 @@ int eea(
   double L,
   _Complex double* u0,
   _Complex double* g,
+  bool irreversible,
   unsigned int print_freq,
   unsigned int starting_time,
   unsigned int nthreads,
@@ -115,10 +117,13 @@ int eea(
 
   // iterate
   for(t=starting_time;t<starting_time+nsteps;t++){
-    ins_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
+    ns_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
+
+    // convolution term
+    ns_T(u,K1,K2,N1,N2,fft1,fft2,ifft);
 
     energy=compute_energy(u, K1, K2);
-    alpha=compute_alpha(u, K1, K2, g);
+    alpha=compute_alpha(u, K1, K2, N1, N2, g, L, ifft.fft);
     enstrophy=compute_enstrophy(u, K1, K2, L);
     
     // running average
@@ -153,6 +158,7 @@ int quiet(
   double L,
   _Complex double* u0,
   _Complex double* g,
+  bool irreversible,
   unsigned int nthreads,
   FILE* savefile
 ){
@@ -171,7 +177,7 @@ int quiet(
 
   // iterate
   for(t=0;t<nsteps;t++){
-    ins_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
+    ns_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
   }
 
   // save final entry to savefile
@@ -267,8 +273,8 @@ int copy_u(
   return 0;
 }
 
-// next time step for Irreversible Navier-Stokes equation
-int ins_step(
+// next time step
+int ns_step(
   _Complex double* u,
   int K1,
   int K2,
@@ -283,12 +289,13 @@ int ins_step(
   fft_vect ifft,
   _Complex double* tmp1,
   _Complex double* tmp2,
-  _Complex double* tmp3
+  _Complex double* tmp3,
+  bool irreversible
 ){
   int kx,ky;
 
   // k1
-  ins_rhs(tmp1, u, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft);
+  ns_rhs(tmp1, u, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
   // add to output
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
@@ -303,7 +310,7 @@ int ins_step(
     }
   }
   // k2
-  ins_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft);
+  ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
   // add to output
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
@@ -318,7 +325,7 @@ int ins_step(
     }
   }
   // k3
-  ins_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft);
+  ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
   // add to output
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
@@ -333,7 +340,7 @@ int ins_step(
     }
   }
   // k4
-  ins_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft);
+  ns_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
   // add to output
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
@@ -344,8 +351,8 @@ int ins_step(
   return(0);
 }
 
-// right side of Irreversible Navier-Stokes equation
-int ins_rhs(
+// right side of Irreversible/Reversible Navier-Stokes equation
+int ns_rhs(
   _Complex double* out,
   _Complex double* u,
   int K1,
@@ -357,14 +364,21 @@ int ins_rhs(
   _Complex double* g,
   fft_vect fft1,
   fft_vect fft2,
-  fft_vect ifft
+  fft_vect ifft,
+  bool irreversible
 ){
   int kx,ky;
   int i;
+  double alpha;
 
   // compute convolution term
   ns_T(u,K1,K2,N1,N2,fft1,fft2,ifft);
 
+  if (irreversible) {
+    alpha=nu;
+  } else {
+    alpha=compute_alpha(u,K1,K2,N1,N2,g,L,ifft.fft);
+  }
 
   /*
   // compare convolution term (store result in fft1.fft)
@@ -376,14 +390,13 @@ int ins_rhs(
   printf("% .15e\n",sqrt(cmp));
   */
 
-
   for(i=0; i<(2*K1+1)*(2*K2+1); i++){
     out[i]=0;
   }
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
       if(kx!=0 || ky!=0){
-        out[klookup(kx,ky,2*K1+1,2*K2+1)]=-4*M_PI*M_PI/L/L*nu*(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]+g[klookup(kx,ky,2*K1+1,2*K2+1)]+4*M_PI*M_PI/L/L/sqrt(kx*kx+ky*ky)*ifft.fft[klookup(kx,ky,N1,N2)];
+        out[klookup(kx,ky,2*K1+1,2*K2+1)]=-4*M_PI*M_PI/L/L*alpha*(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]+g[klookup(kx,ky,2*K1+1,2*K2+1)]+4*M_PI*M_PI/L/L/sqrt(kx*kx+ky*ky)*ifft.fft[klookup(kx,ky,N1,N2)];
       }
     }
   }
@@ -519,6 +532,36 @@ double compute_alpha(
   _Complex double* u,
   int K1,
   int K2,
+  int N1,
+  int N2,
+  _Complex double* g,
+  double L,
+  _Complex double* T
+){
+  _Complex double num=0;
+  double denom=0;
+  int kx,ky;
+
+  num=0.;
+  denom=0.;
+
+  for(kx=-K1;kx<=K1;kx++){
+    for(ky=-K2;ky<=K2;ky++){
+      num+=(L*L/4/M_PI/M_PI*(kx*kx+ky*ky)*g[klookup(kx,ky,2*K1+1,2*K2+1)]+sqrt(kx*kx+ky*ky)*T[klookup(kx,ky,N1,N2)])*conj(u[klookup(kx,ky,2*K1+1,2*K2+1)]);
+      denom+=__real__ (kx*kx+ky*ky)*(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]*conj(u[klookup(kx,ky,2*K1+1,2*K2+1)]);
+    }
+  }
+
+  return __real__ num/denom;
+}
+
+
+/*
+// compute alpha
+double compute_alpha(
+  _Complex double* u,
+  int K1,
+  int K2,
   _Complex double* g
 ){
   _Complex double num=0;
@@ -533,7 +576,7 @@ double compute_alpha(
   }
 
   return __real__ (num/denom);
-}
+}*/
 
 // compute energy
 double compute_energy(