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 +++++++++++++++++++++++++++++++++++++++ src/init.h | 10 +++++++ src/main.c | 84 ++++++++++++++++++++++++++++++++++++++++++++++++----- src/navier-stokes.c | 70 ++++++++++++-------------------------------- src/navier-stokes.h | 12 ++++---- 5 files changed, 172 insertions(+), 66 deletions(-) create mode 100644 src/init.c create mode 100644 src/init.h (limited to 'src') 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; +} 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 #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, ¶m_str, &driving, &command, &nthreads); + ret=read_args(argc, argv, ¶m_str, &driving, &init, &command, &nthreads); if(ret<0){ return(-1); } + // read params ret=read_params(param_str, ¶meters); 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] \n\n"); + fprintf(stderr, "usage:\n nstrophy [-t nthreads] [-p parameters] [-g driving_force] [-i initial_condition] \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;t0){ - 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); -- cgit v1.2.3-54-g00ecf