From 01f47ace6756c28deb9ea0daaee3904ffa5ce9e0 Mon Sep 17 00:00:00 2001 From: Ian Jauslin Date: Thu, 11 Jan 2018 22:48:14 +0000 Subject: Initial commit --- src/main.c | 282 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 282 insertions(+) create mode 100644 src/main.c (limited to 'src/main.c') diff --git a/src/main.c b/src/main.c new file mode 100644 index 0000000..9420e90 --- /dev/null +++ b/src/main.c @@ -0,0 +1,282 @@ +#define VERSION "0.0" + +#include +#include +#include +#include +#include +#include "navier-stokes.h" + +// usage message +int print_usage(); +// read command line arguments +int read_args(int argc, const char* argv[], ns_params* params, unsigned int* nsteps, unsigned int* computation_nr); + +// compute enstrophy as a function of time in the I-NS equation +int enstrophy(ns_params params, unsigned int Nsteps); + + +#define COMPUTATION_ENSTROPHY 1 +int main (int argc, const char* argv[]){ + ns_params params; + unsigned int nsteps; + int ret; + unsigned int computation_nr; + + // default computation: phase diagram + computation_nr=COMPUTATION_ENSTROPHY; + + // read command line arguments + ret=read_args(argc, argv, ¶ms, &nsteps, &computation_nr); + if(ret<0){ + return(-1); + } + if(ret>0){ + return(0); + } + + // enstrophy + if(computation_nr==COMPUTATION_ENSTROPHY){ + enstrophy(params, nsteps); + } + + return(0); +} + +// usage message +int print_usage(){ + fprintf(stderr, "usage:\n nstrophy enstrophy [-h timestep] [-K modes] [-v] [-N nsteps]\n\n nstrophy -V [-v]\n\n"); + return(0); +} + +// read command line arguments +#define CP_FLAG_TIMESTEP 1 +#define CP_FLAG_NSTEPS 2 +#define CP_FLAG_MODES 2 +#define CP_FLAG_NU 3 +int read_args(int argc, const char* argv[], ns_params* params, unsigned int* nsteps, unsigned int* computation_nr){ + int i; + int ret; + // temporary int + int tmp_int; + // temporary unsigned int + unsigned int tmp_uint; + // temporary double + double tmp_double; + // pointers + char* ptr; + // flag that indicates what argument is being read + int flag=0; + // print version and exit + char Vflag=0; + + // defaults + /* + params->h=6.103515625e-05; + params->K=16; + *nsteps=16777216; + params->nu=4.9632717887631524e-05; + */ + params->h=1e-5; + params->K=16; + *nsteps=10000000; + params->nu=1e-4; + + // loop over arguments + for(i=1;ih=tmp_double; + flag=0; + } + // nsteps + else if(flag==CP_FLAG_NSTEPS){ + ret=sscanf(argv[i],"%u",&tmp_uint); + if(ret!=1){ + fprintf(stderr, "error: '-N' should be followed by an unsigned int\n got '%s'\n",argv[i]); + return(-1); + } + *nsteps=tmp_uint; + flag=0; + } + // nsteps + else if(flag==CP_FLAG_MODES){ + ret=sscanf(argv[i],"%d",&tmp_int); + if(ret!=1){ + fprintf(stderr, "error: '-K' should be followed by an int\n got '%s'\n",argv[i]); + return(-1); + } + params->K=tmp_uint; + flag=0; + } + // friction + else if(flag==CP_FLAG_TIMESTEP){ + ret=sscanf(argv[i],"%lf",&tmp_double); + if(ret!=1){ + fprintf(stderr, "error: '-n' should be followed by a double\n got '%s'\n",argv[i]); + return(-1); + } + params->nu=tmp_double; + flag=0; + } + // computation to run + else{ + if(strcmp(argv[i], "enstrophy")==0){ + *computation_nr=COMPUTATION_ENSTROPHY; + } + else{ + fprintf(stderr, "error: unrecognized computation: '%s'\n",argv[i]); + print_usage(); + return(-1); + } + flag=0; + } + } + + // print version and exit + if(Vflag==1){ + printf("nstrophy " VERSION "\n"); + return(1); + } + + return(0); +} + +// compute enstrophy as a function of time in the I-NS equation +int enstrophy(ns_params params, unsigned int Nsteps){ + _Complex double* u; + _Complex double* tmp1; + _Complex double* tmp2; + _Complex double* tmp3; + _Complex double alpha; + _Complex double avg; + unsigned int t; + int kx,ky; + fft_vects fft_vects; + + // sizes + params.S=2*params.K+1; + params.N=4*params.K+1; + + // velocity field + u=calloc(sizeof(_Complex double),params.S*params.S); + params.g=calloc(sizeof(_Complex double),params.S*params.S); + // allocate tmp vectors for computation + tmp1=calloc(sizeof(_Complex double),params.S*params.S); + tmp2=calloc(sizeof(_Complex double),params.S*params.S); + tmp3=calloc(sizeof(_Complex double),params.S*params.S); + + // initial value + for(kx=-params.K;kx<=params.K;kx++){ + for(ky=-params.K;ky<=params.K;ky++){ + //u[KLOOKUP(kx,ky,params.S)]=kx*kx*ky*ky*exp(-(kx*kx+ky*ky)); + if((kx==1 && ky==0) || (kx==-1 && ky==0)){ + u[KLOOKUP(kx,ky,params.S)]=1; + } + else{ + u[KLOOKUP(kx,ky,params.S)]=0; + } + } + } + + // driving force + for(kx=-params.K;kx<=params.K;kx++){ + for(ky=-params.K;ky<=params.K;ky++){ + //params.g[KLOOKUP(kx,ky,params.S)]=sqrt(kx*kx*ky*ky)*exp(-(kx*kx+ky*ky)); + if((kx==2 && ky==-1) || (kx==-2 && ky==1)){ + params.g[KLOOKUP(kx,ky,params.S)]=1; + } + else{ + params.g[KLOOKUP(kx,ky,params.S)]=0; + } + } + } + + + // prepare vectors for fft + fft_vects.fft1=fftw_malloc(sizeof(fftw_complex)*params.N*params.N); + fft_vects.fft1_plan=fftw_plan_dft_2d((int)params.N,(int)params.N, fft_vects.fft1, fft_vects.fft1, FFTW_FORWARD, FFTW_MEASURE); + fft_vects.fft2=fftw_malloc(sizeof(fftw_complex)*params.N*params.N); + fft_vects.fft2_plan=fftw_plan_dft_2d((int)params.N,(int)params.N, fft_vects.fft2, fft_vects.fft2, FFTW_FORWARD, FFTW_MEASURE); + fft_vects.invfft=fftw_malloc(sizeof(fftw_complex)*params.N*params.N); + fft_vects.invfft_plan=fftw_plan_dft_2d((int)params.N,(int)params.N, fft_vects.invfft, fft_vects.invfft, FFTW_BACKWARD, FFTW_MEASURE); + + // init running average + avg=0; + + // iterate + for(t=0;t0){ + avg=avg-(avg-alpha)/t; + } + + if(t>0 && t%1000==0){ + fprintf(stderr,"%8d % .8e % .8e % .8e % .8e\n",t, __real__ avg, __imag__ avg, __real__ alpha, __imag__ alpha); + printf("%8d % .8e % .8e % .8e % .8e\n",t, __real__ avg, __imag__ avg, __real__ alpha, __imag__ alpha); + } + } + + // free memory + fftw_destroy_plan(fft_vects.fft1_plan); + fftw_destroy_plan(fft_vects.fft2_plan); + fftw_destroy_plan(fft_vects.invfft_plan); + fftw_free(fft_vects.fft1); + fftw_free(fft_vects.fft2); + fftw_free(fft_vects.invfft); + + free(tmp3); + free(tmp2); + free(tmp1); + free(params.g); + free(u); + + return(0); +} -- cgit v1.2.3-54-g00ecf