choose initial condition on cli

1 files changed, 62 insertions, 0 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;
+}