From d4254c6b8e6d4a94a5448e771d8a0620f266cf05 Mon Sep 17 00:00:00 2001 From: Ian Jauslin Date: Thu, 26 May 2022 15:05:30 -0400 Subject: choose initial condition on cli --- src/init.c | 62 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 62 insertions(+) create mode 100644 src/init.c (limited to 'src/init.c') 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 +#include + +// 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; +} -- cgit v1.2.3-54-g00ecf