choose initial condition on cli

5 files changed, 172 insertions, 66 deletions

diff --git a/src/init.c b/src/init.c
new file mode 100644
index 0000000..9711db3
--- /dev/null
+++ b/src/init.c
@@ -0,0 +1,62 @@
+#include "init.h"
+#include "navier-stokes.h"
+#include <stdlib.h>
+#include <math.h>
+
+// random initial condition
+int init_random (
+  _Complex double* u0,
+  int K1,
+  int K2,
+  int seed
+){
+  int kx,ky;
+  double rescale;
+  double x,y;
+
+  srand(seed);
+
+  // random init (set half, then the other half are the conjugates)
+  for(kx=0;kx<=K1;kx++){
+    for(ky=-K2;ky<=K2;ky++){
+      if (kx!=0 || ky>0){
+	x=-0.5+((double) rand())/RAND_MAX;
+	y=-0.5+((double) rand())/RAND_MAX;
+	u0[klookup(kx,ky,2*K1+1,2*K2+1)]=x+y*I;
+	u0[klookup(-kx,-ky,2*K1+1,2*K2+1)]=conj(u0[klookup(kx,ky,2*K1+1,2*K2+1)]);
+      }
+    }
+  }
+
+  // rescale to match with Gallavotti's initialization
+  rescale=0;
+  for(kx=-K1;kx<=K1;kx++){
+    for(ky=-K2;ky<=K2;ky++){
+      rescale=rescale+((__real__ u0[klookup(kx,ky,2*K1+1,2*K2+1)])*(__real__ u0[klookup(kx,ky,2*K1+1,2*K2+1)])+(__imag__ u0[klookup(kx,ky,2*K1+1,2*K2+1)])*(__imag__ u0[klookup(kx,ky,2*K1+1,2*K2+1)]))*(kx*kx+ky*ky);
+    }
+  }
+  for(kx=-K1;kx<=K1;kx++){
+    for(ky=-K2;ky<=K2;ky++){
+      u0[klookup(kx,ky,2*K1+1,2*K2+1)]=u0[klookup(kx,ky,2*K1+1,2*K2+1)]*sqrt(1.54511597324389e+02/rescale);
+    }
+  }
+
+  return 0;
+}
+
+// Gaussian initial condition
+int init_gaussian (
+  _Complex double* u0,
+  int K1,
+  int K2
+){
+  int kx,ky;
+
+  for(kx=-K1;kx<=K1;kx++){
+    for(ky=-K2;ky<=K2;ky++){
+      u0[klookup(kx,ky,2*K1+1,2*K2+1)]=(kx*kx+ky*ky)*exp(-(kx*kx+ky*ky));
+    }
+  }
+
+  return 0;
+}
diff --git a/src/init.h b/src/init.h
new file mode 100644
index 0000000..8631b38
--- /dev/null
+++ b/src/init.h
@@ -0,0 +1,10 @@
+#ifndef INIT_H
+#define INIT_H
+
+// random initial condition
+int init_random(_Complex double* u0, int K1, int K2, int seed);
+
+// Gaussian initial condition
+int init_gaussian(_Complex double* u0, int K1, int K2);
+
+#endif
diff --git a/src/main.c b/src/main.c
index 531b9b5..27ed314 100644
--- a/src/main.c
+++ b/src/main.c
@@ -8,6 +8,7 @@
 #include <stdbool.h>
 #include "navier-stokes.h"
 #include "driving.h"
+#include "init.h"
 
 // structure to store parameters, to make it easier to pass parameters to CLI functions
 typedef struct nstrophy_parameters {
@@ -20,16 +21,20 @@ typedef struct nstrophy_parameters {
   double delta;
   double L;
   unsigned int print_freq;
+  int seed;
 } nstrophy_parameters;
 
 // usage message
 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_args(int argc, const char* argv[], char** params, unsigned int* driving_force, unsigned int* command, unsigned int* init, unsigned int* nthreads);
 int read_params(char* param_str, nstrophy_parameters* parameters);
 int set_parameter(char* lhs, char* rhs, nstrophy_parameters* parameters, bool* setN1, bool* setN2);
 
+// set initial condition
+_Complex double* set_init(unsigned int init, nstrophy_parameters parameters);
+
 #define COMMAND_UK 1
 #define COMMAND_ENSTROPHY 2
 #define COMMAND_QUIET 3
@@ -38,6 +43,9 @@ int set_parameter(char* lhs, char* rhs, nstrophy_parameters* parameters, bool* s
 #define DRIVING_ZERO 1
 #define DRIVING_TEST 2
 
+#define INIT_RANDOM 1
+#define INIT_GAUSSIAN 2
+
 
 int main (
   int argc,
@@ -47,22 +55,26 @@ int main (
   nstrophy_parameters parameters;
   _Complex double (*g)(int,int);
   int ret;
-  unsigned int driving,command;
+  unsigned int driving,command,init;
   unsigned int nthreads=1;
+  _Complex double* u0;
 
   command=0;
   driving=0;
+  init=0;
 
   // read command line arguments
-  ret=read_args(argc, argv, &param_str, &driving, &command, &nthreads);
+  ret=read_args(argc, argv, &param_str, &driving, &init, &command, &nthreads);
   if(ret<0){
     return(-1);
   }
+
   // read params
   ret=read_params(param_str, &parameters);
   if(ret<0){
     return(-1);
   }
+
   // set driving force
   switch(driving){
   case DRIVING_ZERO:
@@ -77,30 +89,35 @@ int main (
     break;
   }
 
+  // set initial condition
+  u0=set_init(init, parameters);
+
   // run command
   if (command==COMMAND_UK){
-    uk(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, g, parameters.print_freq, nthreads);
+    uk(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.print_freq, nthreads);
   }
   else if (command==COMMAND_ENERGY){
-    energy(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, g, parameters.print_freq, nthreads);
+    energy(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.print_freq, nthreads);
   }
   else if(command==COMMAND_ENSTROPHY){
-    enstrophy(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, g, parameters.print_freq, nthreads);
+    enstrophy(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.print_freq, nthreads);
   }
   else if(command==COMMAND_QUIET){
-    quiet(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, g, nthreads);
+    quiet(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, nthreads);
   }
   else if(command==0){
     fprintf(stderr, "error: no command specified\n");
     print_usage();
   }
 
+  free(u0);
+
   return(0);
 }
 
 // usage message
 int print_usage(){
-  fprintf(stderr, "usage:\n       nstrophy [-p parameters] [-g driving_force] <command>\n\n");
+  fprintf(stderr, "usage:\n       nstrophy [-t nthreads] [-p parameters] [-g driving_force] [-i initial_condition] <command>\n\n");
   return(0);
 }
 
@@ -108,11 +125,13 @@ int print_usage(){
 #define CP_FLAG_PARAMS 1
 #define CP_FLAG_DRIVING 2
 #define CP_FLAG_NTHREADS 3
+#define CP_FLAG_INIT 4
 int read_args(
   int argc,
   const char* argv[],
   char** params,
   unsigned int* driving_force,
+  unsigned int* init,
   unsigned int* command,
   unsigned int* nthreads
 ){
@@ -138,6 +157,9 @@ int read_args(
 	case 't':
 	  flag=CP_FLAG_NTHREADS;
 	  break;
+	case 'i':
+	  flag=CP_FLAG_INIT;
+	  break;
 	default:
 	  fprintf(stderr, "unrecognized option '-%c'\n", *ptr);
 	  print_usage();
@@ -174,6 +196,20 @@ int read_args(
       }
       flag=0;
     }
+    // initial condition
+    else if(flag==CP_FLAG_INIT){
+      if (strcmp(argv[i],"random")==0){
+	*init=INIT_RANDOM;
+      }
+      else if (strcmp(argv[i],"gaussian")==0){
+	*init=INIT_GAUSSIAN;
+      }
+      else{
+	fprintf(stderr, "error: unrecognized initial condition '%s'\n",argv[i]);
+	return(-1);
+      }
+      flag=0;
+    }
     // computation to run
     else{
       if(strcmp(argv[i], "uk")==0){
@@ -226,6 +262,7 @@ int read_params(
   parameters->L=2*M_PI;
   parameters->nsteps=10000000;
   parameters->print_freq=1000;
+  parameters->seed=17;
 
   if (param_str!=NULL){
     // init
@@ -389,6 +426,13 @@ int set_parameter(
       return(-1);
     }
   }
+  else if (strcmp(lhs,"random_seed")==0){
+    ret=sscanf(rhs,"%d",&(parameters->seed));
+    if(ret!=1){
+      fprintf(stderr, "error: parameter 'random_seed' should be an integer\n       got '%s'\n",rhs);
+      return(-1);
+    }
+  }
   else{
     fprintf(stderr, "error: unrecognized parameter '%s'\n",lhs);
     return(-1);
@@ -396,3 +440,27 @@ int set_parameter(
 
   return(0);
 }
+
+// set initial condition
+_Complex double* set_init(
+  unsigned int init,
+  nstrophy_parameters parameters
+){
+  _Complex double* u0=calloc(sizeof(_Complex double),(2*parameters.K1+1)*(2*parameters.K2+1));
+
+  switch(init){
+  case INIT_RANDOM:
+    init_random(u0, parameters.K1, parameters.K2, parameters.seed);
+    break;
+
+  case INIT_GAUSSIAN:
+    init_gaussian(u0, parameters.K1, parameters.K2);
+    break;
+
+  default:
+    init_gaussian(u0, parameters.K1, parameters.K2);
+    break;
+  }
+  
+  return u0;
+}
diff --git a/src/navier-stokes.c b/src/navier-stokes.c
index 798c86a..c52b810 100644
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -12,6 +12,7 @@ int uk(
   double nu,
   double delta,
   double L,
+  _Complex double* u0,
   _Complex double (*g)(int,int),
   unsigned int print_freq,
   unsigned int nthreads
@@ -27,7 +28,8 @@ int uk(
   int kx,ky;
 
   ns_init_tmps(&u, &tmp1, &tmp2, &tmp3, &fft1, &fft2, &ifft, K1, K2, N1, N2, nthreads);
-  ns_init_u(u, K1, K2);
+  // copy initial condition
+  copy_u(u, u0, K1, K2);
 
   // print column headers
   printf("#    1:i                    2:t ");
@@ -77,6 +79,7 @@ int energy(
   double nu,
   double delta,
   double L,
+  _Complex double* u0,
   _Complex double (*g)(int,int),
   unsigned int print_freq,
   unsigned int nthreads
@@ -93,7 +96,8 @@ int energy(
   double energy;
 
   ns_init_tmps(&u, &tmp1, &tmp2, &tmp3, &fft1, &fft2, &ifft, K1, K2, N1, N2, nthreads);
-  ns_init_u(u, K1, K2);
+  // copy initial condition
+  copy_u(u, u0, K1, K2);
 
   // iterate
   for(t=0;t<nsteps;t++){
@@ -126,6 +130,7 @@ int enstrophy(
   double nu,
   double delta,
   double L,
+  _Complex double* u0,
   _Complex double (*g)(int,int),
   unsigned int print_freq,
   unsigned int nthreads
@@ -142,7 +147,8 @@ int enstrophy(
   fft_vect ifft;
 
   ns_init_tmps(&u, &tmp1, &tmp2, &tmp3, &fft1, &fft2, &ifft, K1, K2, N1, N2, nthreads);
-  ns_init_u(u, K1, K2);
+  // copy initial condition
+  copy_u(u, u0, K1, K2);
 
 
   // init running average
@@ -178,6 +184,7 @@ int quiet(
   double nu,
   double delta,
   double L,
+  _Complex double* u0,
   _Complex double (*g)(int,int),
   unsigned int nthreads
 ){
@@ -191,7 +198,8 @@ int quiet(
   fft_vect ifft;
 
   ns_init_tmps(&u, &tmp1, &tmp2, &tmp3, &fft1, &fft2, &ifft, K1, K2, N1, N2, nthreads);
-  ns_init_u(u, K1, K2);
+  // copy initial condition
+  copy_u(u, u0, K1, K2);
 
   // iterate
   for(t=0;t<nsteps;t++){
@@ -272,59 +280,17 @@ int ns_free_tmps(
 
 
 
-// initial value
-int ns_init_u(
+// copy u0 to u
+int copy_u(
   _Complex double* u,
+  _Complex double* u0,
   int K1,
   int K2
 ){
-  int kx,ky;
-
-  /*
-  srand(17);
-
-  // random init (set half, then the other half are the conjugates)
-  for(kx=0;kx<=K1;kx++){
-    for(ky=-K2;ky<=K2;ky++){
-      if (kx!=0 || ky>0){
-	double x=-0.5+((double) rand())/RAND_MAX;
-	double y=-0.5+((double) rand())/RAND_MAX;
-	u[klookup(kx,ky,2*K1+1,2*K2+1)]=x+y*I;
-	u[klookup(-kx,-ky,2*K1+1,2*K2+1)]=conj(u[klookup(kx,ky,2*K1+1,2*K2+1)]);
-      }
-    }
-  }
-
-  // rescale to match with Gallavotti's initialization
-  double rescale;
-  rescale=0;
-  for(kx=-K1;kx<=K1;kx++){
-    for(ky=-K2;ky<=K2;ky++){
-      rescale=rescale+((__real__ u[klookup(kx,ky,2*K1+1,2*K2+1)])*(__real__ u[klookup(kx,ky,2*K1+1,2*K2+1)])+(__imag__ u[klookup(kx,ky,2*K1+1,2*K2+1)])*(__imag__ u[klookup(kx,ky,2*K1+1,2*K2+1)]))*(kx*kx+ky*ky);
-    }
-  }
-  for(kx=-K1;kx<=K1;kx++){
-    for(ky=-K2;ky<=K2;ky++){
-      u[klookup(kx,ky,2*K1+1,2*K2+1)]=u[klookup(kx,ky,2*K1+1,2*K2+1)]*sqrt(1.54511597324389e+02/rescale);
-    }
-  }
-  */
-
-
-  /*
-  // constant init
-  for(kx=-K1;kx<=K1;kx++){
-    for(ky=-K2;ky<=K2;ky++){
-      u[klookup(kx,ky,2*K1+1,2*K2+1)]=1.;
-    }
-  }
-  */
+  int i;
 
-  // gaussian init
-  for(kx=-K1;kx<=K1;kx++){
-    for(ky=-K2;ky<=K2;ky++){
-      u[klookup(kx,ky,2*K1+1,2*K2+1)]=(kx*kx+ky*ky)*exp(-(kx*kx+ky*ky));
-    }
+  for(i=0;i<(2*K1+1)*(2*K2+1);i++){
+    u[i]=u0[i];
   }
 
   return 0;
diff --git a/src/navier-stokes.h b/src/navier-stokes.h
index 714df7d..e23b457 100644
--- a/src/navier-stokes.h
+++ b/src/navier-stokes.h
@@ -13,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, double L, _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* u0, _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, double L, _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* u0, _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, double L, _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* u0, _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, double L, _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* u0, _Complex double (*g)(int,int), unsigned int nthreads);
 
 
 // initialize vectors for computation
@@ -30,8 +30,8 @@ int  ns_init_tmps( _Complex double **u, _Complex double ** tmp1, _Complex double
 // release vectors
 int ns_free_tmps( _Complex double* u, _Complex double* tmp1, _Complex double *tmp2,_Complex double *tmp3, fft_vect fft1, fft_vect fft2, fft_vect ifft);
 
-// initial value
-int ns_init_u( _Complex double* u, int K1, int K2);
+// 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)(int,int), fft_vect fft1, fft_vect fft2,fft_vect ifft, _Complex double* tmp1, _Complex double *tmp2, _Complex double *tmp3);