add L parameter

3 files changed, 43 insertions, 28 deletions

diff --git a/src/main.c b/src/main.c
index 27bfb7a..3d8867e 100644
--- a/src/main.c
+++ b/src/main.c
@@ -14,8 +14,8 @@ int print_usage();
 
 // read command line arguments
 int read_args(int argc, const char* argv[], char** params, unsigned int* driving_force, unsigned int* command, unsigned int* nthreads);
-int read_params(char* params, int* K1, int* K2, int* N1, int* N2, unsigned int* nsteps, double* nu, double* delta, unsigned int* print_freq);
-int set_parameter(char* lhs, char* rhs, int* K1, int* K2, int* N1, int* N2, unsigned int* nsteps, double* nu, double* delta, unsigned int* print_freq, bool* setN1, bool* setN2);
+int read_params(char* params, int* K1, int* K2, int* N1, int* N2, unsigned int* nsteps, double* nu, double* delta, double* L, unsigned int* print_freq);
+int set_parameter(char* lhs, char* rhs, int* K1, int* K2, int* N1, int* N2, unsigned int* nsteps, double* nu, double* delta, double* L, unsigned int* print_freq, bool* setN1, bool* setN2);
 
 
 #define COMMAND_UK 1
@@ -35,7 +35,7 @@ int main (
   int K1,K2;
   int N1,N2;
   unsigned int nsteps;
-  double nu,delta;
+  double nu,delta,L;
   _Complex double (*g)(int,int);
   int ret;
   unsigned int driving,command;
@@ -51,7 +51,7 @@ int main (
     return(-1);
   }
   // read params
-  ret=read_params(params, &K1, &K2, &N1, &N2, &nsteps, &nu, &delta, &print_freq);
+  ret=read_params(params, &K1, &K2, &N1, &N2, &nsteps, &nu, &delta, &L, &print_freq);
   if(ret<0){
     return(-1);
   }
@@ -71,16 +71,16 @@ int main (
 
   // run command
   if (command==COMMAND_UK){
-    uk(K1, K2, N1, N2, nsteps, nu, delta, g, print_freq, nthreads);
+    uk(K1, K2, N1, N2, nsteps, nu, delta, L, g, print_freq, nthreads);
   }
   else if (command==COMMAND_ENERGY){
-    energy(K1, K2, N1, N2, nsteps, nu, delta, g, print_freq, nthreads);
+    energy(K1, K2, N1, N2, nsteps, nu, delta, L, g, print_freq, nthreads);
   }
   else if(command==COMMAND_ENSTROPHY){
-    energy(K1, K2, N1, N2, nsteps, nu, delta, g, print_freq, nthreads);
+    energy(K1, K2, N1, N2, nsteps, nu, delta, L, g, print_freq, nthreads);
   }
   else if(command==COMMAND_QUIET){
-    quiet(K1, K2, N1, N2, nsteps, nu, delta, g, nthreads);
+    quiet(K1, K2, N1, N2, nsteps, nu, delta, L, g, nthreads);
   }
   else if(command==0){
     fprintf(stderr, "error: no command specified\n");
@@ -201,6 +201,7 @@ int read_params(
   unsigned int* nsteps,
   double* nu,
   double* delta,
+  double* L,
   unsigned int* print_freq
 ){
   int ret;
@@ -222,6 +223,7 @@ int read_params(
   *delta=0.0001220703125;
   //nu=2^-11
   *nu=0.00048828125;
+  *L=2*M_PI;
   *nsteps=10000000;
   *print_freq=1000;
 
@@ -244,7 +246,7 @@ int read_params(
 	break;
       case ';':
 	//set parameter
-	ret=set_parameter(buffer_lhs,buffer_rhs,K1,K2,N1,N2,nsteps,nu,delta,print_freq,&setN1,&setN2);
+	ret=set_parameter(buffer_lhs,buffer_rhs,K1,K2,N1,N2,nsteps,nu,delta,L,print_freq,&setN1,&setN2);
 	if(ret<0){
 	  return ret;
 	}
@@ -271,7 +273,7 @@ int read_params(
 
     // set last param
     if (*params!='\0'){
-      ret=set_parameter(buffer_lhs,buffer_rhs,K1,K2,N1,N2,nsteps,nu,delta,print_freq,&setN1,&setN2);
+      ret=set_parameter(buffer_lhs,buffer_rhs,K1,K2,N1,N2,nsteps,nu,delta,L,print_freq,&setN1,&setN2);
       if(ret<0){
 	return ret;
       }
@@ -305,6 +307,7 @@ int set_parameter(
   unsigned int* nsteps,
   double* nu,
   double* delta,
+  double* L,
   unsigned int* print_freq,
   bool* setN1,
   bool* setN2
@@ -380,6 +383,13 @@ int set_parameter(
       return(-1);
     }
   }
+  else if (strcmp(lhs,"L")==0){
+    ret=sscanf(rhs,"%lf",L);
+    if(ret!=1){
+      fprintf(stderr, "error: parameter 'L' should be a double\n       got '%s'\n",rhs);
+      return(-1);
+    }
+  }
   else if (strcmp(lhs,"print_freq")==0){
     ret=sscanf(rhs,"%u",print_freq);
     if(ret!=1){
diff --git a/src/navier-stokes.c b/src/navier-stokes.c
index 99e41ab..89bfbbd 100644
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -2,8 +2,6 @@
 #include <math.h>
 #include <stdlib.h>
 
-#define M_PI 3.14159265358979323846
-
 // compute solution as a function of time
 int uk(
   int K1,
@@ -13,6 +11,7 @@ int uk(
   unsigned int nsteps,
   double nu,
   double delta,
+  double L,
   _Complex double (*g)(int,int),
   unsigned int print_freq,
   unsigned int nthreads
@@ -44,7 +43,7 @@ int uk(
 
   // iterate
   for(t=0;t<nsteps;t++){
-    ins_step(u, K1, K2, N1, N2, nu, delta, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
+    ins_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
     
     if(t%print_freq==0){
       fprintf(stderr,"%d % .8e ",t,t*delta);
@@ -77,6 +76,7 @@ int energy(
   unsigned int nsteps,
   double nu,
   double delta,
+  double L,
   _Complex double (*g)(int,int),
   unsigned int print_freq,
   unsigned int nthreads
@@ -97,7 +97,7 @@ int energy(
 
   // iterate
   for(t=0;t<nsteps;t++){
-    ins_step(u, K1, K2, N1, N2, nu, delta, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
+    ins_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
     
     if(t%print_freq==0){
 
@@ -126,6 +126,7 @@ int enstrophy(
   unsigned int nsteps,
   double nu,
   double delta,
+  double L,
   _Complex double (*g)(int,int),
   unsigned int print_freq,
   unsigned int nthreads
@@ -150,7 +151,7 @@ int enstrophy(
 
   // iterate
   for(t=0;t<nsteps;t++){
-    ins_step(u, K1, K2, N1, N2, nu, delta, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
+    ins_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
     alpha=compute_alpha(u, K1, K2, g);
     
     // running average
@@ -177,6 +178,7 @@ int quiet(
   unsigned int nsteps,
   double nu,
   double delta,
+  double L,
   _Complex double (*g)(int,int),
   unsigned int nthreads
 ){
@@ -194,7 +196,7 @@ int quiet(
 
   // iterate
   for(t=0;t<nsteps;t++){
-    ins_step(u, K1, K2, N1, N2, nu, delta, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
+    ins_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
   }
 
   ns_free_tmps(u, tmp1, tmp2, tmp3, fft1, fft2, ifft);
@@ -342,6 +344,7 @@ int ins_step(
   int N2,
   double nu,
   double delta,
+  double L,
   _Complex double (*g)(int,int),
   fft_vect fft1,
   fft_vect fft2,
@@ -353,7 +356,7 @@ int ins_step(
   int kx,ky;
 
   // k1
-  ins_rhs(tmp1, u, K1, K2, N1, N2, nu, g, fft1, fft2, ifft);
+  ins_rhs(tmp1, u, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft);
   // add to output
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
@@ -368,7 +371,7 @@ int ins_step(
     }
   }
   // k2
-  ins_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, g, fft1, fft2, ifft);
+  ins_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft);
   // add to output
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
@@ -383,7 +386,7 @@ int ins_step(
     }
   }
   // k3
-  ins_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, g, fft1, fft2, ifft);
+  ins_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft);
   // add to output
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
@@ -398,7 +401,7 @@ int ins_step(
     }
   }
   // k4
-  ins_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, g, fft1, fft2, ifft);
+  ins_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft);
   // add to output
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
@@ -418,6 +421,7 @@ int ins_rhs(
   int N1,
   int N2,
   double nu,
+  double L,
   _Complex double (*g)(int,int),
   fft_vect fft1,
   fft_vect fft2,
@@ -447,7 +451,7 @@ int ins_rhs(
   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*nu*(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]+(*g)(kx,ky)+4*M_PI*M_PI/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*nu*(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]+(*g)(kx,ky)+4*M_PI*M_PI/L/L/sqrt(kx*kx+ky*ky)*ifft.fft[klookup(kx,ky,N1,N2)];
       }
     }
   }
diff --git a/src/navier-stokes.h b/src/navier-stokes.h
index 421342a..714df7d 100644
--- a/src/navier-stokes.h
+++ b/src/navier-stokes.h
@@ -4,6 +4,7 @@
 #include <complex.h>
 #include <fftw3.h>
 
+#define M_PI 3.14159265358979323846
 
 // arrays on which the ffts are performed
 typedef struct fft_vects {
@@ -12,16 +13,16 @@ 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, _Complex double (*g)(int,int), unsigned int print_freq, unsigned int nthreads);
+int uk( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double (*g)(int,int), unsigned int print_freq, unsigned int nthreads);
 
 // compute the energy as a function of time
-int energy( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, _Complex double (*g)(int,int), unsigned int print_freq, unsigned int nthreads);
+int energy( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double (*g)(int,int), unsigned int print_freq, unsigned int nthreads);
 
 // compute enstrophy
-int enstrophy( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, _Complex double (*g)(int,int), unsigned int print_freq, unsigned int nthreads);
+int enstrophy( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double (*g)(int,int), unsigned int print_freq, unsigned int nthreads);
 
 // 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, _Complex double (*g)(int,int), unsigned int nthreads);
+int quiet( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double (*g)(int,int), unsigned int nthreads);
 
 
 // initialize vectors for computation
@@ -33,15 +34,15 @@ int ns_free_tmps( _Complex double* u, _Complex double* tmp1, _Complex double *tm
 int ns_init_u( _Complex double* u, 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, _Complex double (*g)(int,int), fft_vect fft1, fft_vect fft2,fft_vect ifft, _Complex double* tmp1, _Complex double *tmp2, _Complex double *tmp3);
+int ins_step( _Complex double* u, int K1, int K2, int N1, int N2, double nu, double delta, double L, _Complex double (*g)(int,int), fft_vect fft1, fft_vect fft2,fft_vect ifft, _Complex double* tmp1, _Complex double *tmp2, _Complex double *tmp3);
 
 // 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, _Complex double (*g)(int,int), fft_vect fft1, fft_vect fft2, fft_vect ifft);
+int ins_rhs( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2, double nu, double L, _Complex double (*g)(int,int), fft_vect fft1, fft_vect fft2, fft_vect ifft);
 
 // 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);
 
-// convolution term in right side of equation
+// convolution term in right side of equation (computed without fft)
 int ns_T_nofft( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2);
 
 // compute alpha