choose initial condition on cli

1 files changed, 18 insertions, 52 deletions

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;