Reversible equation

3 files changed, 87 insertions, 29 deletions

diff --git a/src/main.c b/src/main.c
index 853f172..f60a38e 100644
--- a/src/main.c
+++ b/src/main.c
@@ -14,6 +14,7 @@
 
 // structure to store parameters, to make it easier to pass parameters to CLI functions
 typedef struct nstrophy_parameters {
+  bool irreversible;
   int K1;
   int K2;
   int N1;
@@ -120,13 +121,13 @@ int main (
 
   // run command
   if (command==COMMAND_UK){
-    uk(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.print_freq, parameters.starting_time, nthreads, savefile);
+    uk(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.irreversible, parameters.print_freq, parameters.starting_time, nthreads, savefile);
   }
   else if(command==COMMAND_EEA){
-    eea(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.print_freq, parameters.starting_time, nthreads, savefile);
+    eea(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.irreversible, parameters.print_freq, parameters.starting_time, nthreads, savefile);
   }
   else if(command==COMMAND_QUIET){
-    quiet(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, nthreads, savefile);
+    quiet(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.irreversible, nthreads, savefile);
   }
   else if(command==0){
     fprintf(stderr, "error: no command specified\n");
@@ -337,6 +338,7 @@ int read_params(
   bool lhs=true;
 
   // defaults
+  parameters->irreversible=true;
   parameters->K1=16;
   parameters->K2=parameters->K1;
   //delta=2^-13
@@ -428,7 +430,19 @@ int set_parameter(
 ){
   int ret;
 
-  if (strcmp(lhs,"K1")==0){
+  if (strcmp(lhs,"equation")==0){
+    if (strcmp(rhs,"irreversible")==0){
+      parameters->irreversible=true;
+    }
+    else if (strcmp(rhs,"reversible")==0){
+      parameters->irreversible=false;
+    }
+    else {
+      fprintf(stderr, "error: 'equation' should be 'irreversible' or 'reversible'\n       got '%s'\n",rhs);
+      return(-1);
+    }
+  }
+  else if (strcmp(lhs,"K1")==0){
     ret=sscanf(rhs,"%d",&(parameters->K1));
     if(ret!=1){
       fprintf(stderr, "error: parameter 'K1' should be an integer\n       got '%s'\n",rhs);
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(
diff --git a/src/navier-stokes.h b/src/navier-stokes.h
index 57399a4..582212a 100644
--- a/src/navier-stokes.h
+++ b/src/navier-stokes.h
@@ -2,6 +2,7 @@
 #define NAVIERSTOKES_H
 
 #include <complex.h>
+#include <stdbool.h>
 #include <fftw3.h>
 
 #define M_PI 3.14159265358979323846
@@ -13,13 +14,13 @@ typedef struct fft_vects {
 } fft_vect;
 
 // compute u_k
-int uk( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, unsigned int print_freq, unsigned int starting_time, unsigned int nthreadsl, FILE* savefile);
+int uk( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int print_freq, unsigned int starting_time, unsigned int nthreadsl, FILE* savefile);
 
 // compute energy, enstrophy and alpha
-int eea( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, unsigned int print_freq, unsigned int starting_time, unsigned int nthreads, FILE* savefile);
+int eea( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int print_freq, unsigned int starting_time, unsigned int nthreads, FILE* savefile);
 
 // compute solution as a function of time, but do not print anything (useful for debugging)
-int quiet( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, unsigned int nthreads, FILE* savefile);
+int quiet( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int nthreads, FILE* savefile);
 
 
 // initialize vectors for computation
@@ -30,11 +31,11 @@ int ns_free_tmps( _Complex double* u, _Complex double* tmp1, _Complex double *tm
 // copy u0 to u
 int copy_u( _Complex double* u, _Complex double* u0, int K1, int K2);
 
-// next time step for Irreversible Navier-Stokes equation
-int ins_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);
+// 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);
 
-// right side of Irreversible Navier-Stokes equation
-int ins_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);
+// 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);
 
 // convolution term in right side of equation
 int ns_T( _Complex double* u, int K1, int K2, int N1, int N2, fft_vect fft1, fft_vect fft2, fft_vect ifft);
@@ -43,7 +44,7 @@ int ns_T( _Complex double* u, int K1, int K2, int N1, int N2, fft_vect fft1, fft
 int ns_T_nofft( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2);
 
 // compute alpha
-double compute_alpha( _Complex double* u, int K1, int K2, _Complex double* g);
+double compute_alpha( _Complex double* u, int K1, int K2, int N1, int N2, _Complex double* g, double L, _Complex double* T);
 // compute energy
 double compute_energy( _Complex double* u, int K1, int K2);
 // compute enstrophy