diff --git a/README.md b/README.md index fd1b9b2..0689759 100644 --- a/README.md +++ b/README.md @@ -125,6 +125,9 @@ should be a `;` sperated list of `key=value` pairs. The possible keys are f_k/|k|^3, `k32` for the normalized L1 norm, `enstrophy` for the enstrophy norm. +* `keep_en_cst` (0 or 1, default 0): impose that the enstrophy is constant at + each step (only really useful for the reversible equation). + # Interrupting and resuming the computation diff --git a/src/main.c b/src/main.c index d3efda8..eedc1e5 100644 --- a/src/main.c +++ b/src/main.c @@ -50,6 +50,7 @@ typedef struct nstrophy_parameters { unsigned int driving; unsigned int init; unsigned int algorithm; + bool keep_en_cst; FILE* initfile; FILE* drivingfile; } nstrophy_parameters; @@ -178,16 +179,16 @@ int main ( // run command if (command==COMMAND_UK){ - uk(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_norm, u0, g, parameters.irreversible, parameters.algorithm, parameters.print_freq, parameters.starting_time, nthreads, savefile); + uk(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_norm, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_en, parameters.algorithm, parameters.print_freq, parameters.starting_time, nthreads, savefile); } else if(command==COMMAND_ENSTROPHY){ // register signal handler to handle aborts signal(SIGINT, sig_handler); signal(SIGTERM, sig_handler); - enstrophy(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_norm, u0, g, parameters.irreversible, parameters.algorithm, parameters.print_freq, parameters.starting_time, nthreads, savefile, utfile, (char*)argv[0], param_str, savefile_str, utfile_str); + enstrophy(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_norm, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_en, parameters.algorithm, parameters.print_freq, parameters.starting_time, nthreads, savefile, utfile, (char*)argv[0], param_str, savefile_str, utfile_str); } else if(command==COMMAND_QUIET){ - quiet(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_norm, parameters.starting_time, u0, g, parameters.irreversible, parameters.algorithm, nthreads, savefile); + quiet(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_norm, parameters.starting_time, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_en, parameters.algorithm, nthreads, savefile); } else if(command==0){ fprintf(stderr, "error: no command specified\n"); @@ -448,6 +449,7 @@ int read_params( parameters->init=INIT_GAUSSIAN; parameters->initfile=NULL; parameters->algorithm=ALGORITHM_RK4; + parameters->keep_en_cst=false; if (param_str!=NULL){ // init @@ -756,6 +758,15 @@ int set_parameter( return(-1); } } + else if (strcmp(lhs,"keep_en_cst")==0){ + int tmp; + ret=sscanf(rhs,"%d",&tmp); + if(ret!=1 || (tmp!=0 && tmp!=1)){ + fprintf(stderr, "error: parameter 'keep_en_cst' should be 0 or 1\n got '%s'\n",rhs); + return(-1); + } + parameters->keep_en_cst=(tmp==1); + } else{ fprintf(stderr, "error: unrecognized parameter '%s'\n",lhs); return(-1); diff --git a/src/navier-stokes.c b/src/navier-stokes.c index 5af1c13..9ab9126 100644 --- a/src/navier-stokes.c +++ b/src/navier-stokes.c @@ -40,6 +40,8 @@ int uk( _Complex double* u0, _Complex double* g, bool irreversible, + bool keep_en_cst, + double target_en, unsigned int algorithm, double print_freq, double starting_time, @@ -88,17 +90,17 @@ int uk( time=starting_time; while(final_time<0 || time0 ? -K2 : 1);ky<=K2;ky++){ + u[klookup_sym(kx,ky,K2)]*=sqrt(target_en/en); + } + } + } + return(0); } @@ -635,7 +653,9 @@ int ns_step_rk2( fft_vect ifft, _Complex double* tmp1, _Complex double* tmp2, - bool irreversible + bool irreversible, + bool keep_en_cst, + double target_en ){ int kx,ky; @@ -657,6 +677,16 @@ int ns_step_rk2( } } + // keep enstrophy constant + if(keep_en_cst){ + double en=compute_enstrophy(u, K1, K2, L); + for(kx=0;kx<=K1;kx++){ + for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){ + u[klookup_sym(kx,ky,K2)]*=sqrt(target_en/en); + } + } + } + return(0); } @@ -688,6 +718,8 @@ int ns_step_rkf45( _Complex double* k6, _Complex double* tmp, bool irreversible, + bool keep_en_cst, + double target_en, // whether to compute k1 or leave it as is bool compute_k1 ){ @@ -805,12 +837,22 @@ int ns_step_rkf45( } // next delta to use in future steps (do not exceed max_delta) *next_delta=fmin(max_delta, (*delta)*pow(relative*tolerance/err,0.2)); + + // keep enstrophy constant + if(keep_en_cst){ + double en=compute_enstrophy(u, K1, K2, L); + for(kx=0;kx<=K1;kx++){ + for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){ + u[klookup_sym(kx,ky,K2)]*=sqrt(target_en/en); + } + } + } } // error too big: repeat with smaller step else{ *delta=factor*(*delta)*pow(relative*tolerance/err,0.2); // do not recompute k1 - ns_step_rkf45(u,tolerance,factor,max_delta,adaptive_norm,K1,K2,N1,N2,nu,delta,next_delta,L,g,fft1,fft2,ifft,k1,k2,k3,k4,k5,k6,tmp,irreversible,false); + ns_step_rkf45(u,tolerance,factor,max_delta,adaptive_norm,K1,K2,N1,N2,nu,delta,next_delta,L,g,fft1,fft2,ifft,k1,k2,k3,k4,k5,k6,tmp,irreversible,keep_en_cst,target_en,false); } return 0; @@ -844,6 +886,8 @@ int ns_step_rkbs32( _Complex double** k4, _Complex double* tmp, bool irreversible, + bool keep_en_cst, + double target_en, // whether to compute k1 bool compute_k1 ){ @@ -944,12 +988,22 @@ int ns_step_rkbs32( tmp=*k1; *k1=*k4; *k4=tmp; + + // keep enstrophy constant + if(keep_en_cst){ + double en=compute_enstrophy(u, K1, K2, L); + for(kx=0;kx<=K1;kx++){ + for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){ + u[klookup_sym(kx,ky,K2)]*=sqrt(target_en/en); + } + } + } } // error too big: repeat with smaller step else{ *delta=factor*(*delta)*pow(relative*tolerance/err,1./3); // this will reuse the same k1 without re-computing it - ns_step_rkbs32(u,tolerance,factor,max_delta,adaptive_norm,K1,K2,N1,N2,nu,delta,next_delta,L,g,fft1,fft2,ifft,k1,k2,k3,k4,tmp,irreversible,false); + ns_step_rkbs32(u,tolerance,factor,max_delta,adaptive_norm,K1,K2,N1,N2,nu,delta,next_delta,L,g,fft1,fft2,ifft,k1,k2,k3,k4,tmp,irreversible,keep_en_cst,target_en,false); } return 0; @@ -984,6 +1038,8 @@ int ns_step_rkdp54( _Complex double* k6, _Complex double* tmp, bool irreversible, + bool keep_en_cst, + double target_en, // whether to compute k1 bool compute_k1 ){ @@ -1109,12 +1165,22 @@ int ns_step_rkdp54( tmp=*k1; *k1=*k2; *k2=tmp; + + // keep enstrophy constant + if(keep_en_cst){ + double en=compute_enstrophy(u, K1, K2, L); + for(kx=0;kx<=K1;kx++){ + for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){ + u[klookup_sym(kx,ky,K2)]*=sqrt(target_en/en); + } + } + } } // error too big: repeat with smaller step else{ *delta=factor*(*delta)*pow(relative*tolerance/err,1./5); // this will reuse the same k1 without re-computing it - ns_step_rkdp54(u,tolerance,factor,max_delta,adaptive_norm,K1,K2,N1,N2,nu,delta,next_delta,L,g,fft1,fft2,ifft,k1,k2,k3,k4,k5,k6,tmp,irreversible,false); + ns_step_rkdp54(u,tolerance,factor,max_delta,adaptive_norm,K1,K2,N1,N2,nu,delta,next_delta,L,g,fft1,fft2,ifft,k1,k2,k3,k4,k5,k6,tmp,irreversible,keep_en_cst,target_en,false); } return 0; diff --git a/src/navier-stokes.h b/src/navier-stokes.h index 639d361..34a61d0 100644 --- a/src/navier-stokes.h +++ b/src/navier-stokes.h @@ -33,13 +33,13 @@ typedef struct fft_vects { } fft_vect; // compute u_k -int uk( int K1, int K2, int N1, int N2, double final_time, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_norm, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int algorithm, double print_freq, double starting_time, unsigned int nthreadsl, FILE* savefile); +int uk( int K1, int K2, int N1, int N2, double final_time, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_norm, _Complex double* u0, _Complex double* g, bool irreversible, bool keep_en_cst, double target_en, unsigned int algorithm, double print_freq, double starting_time, unsigned int nthreadsl, FILE* savefile); // compute enstrophy and alpha -int enstrophy( int K1, int K2, int N1, int N2, double final_time, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_norm, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int algorithm, double print_freq, double starting_time, unsigned int nthreads, FILE* savefile, FILE* utfile, char* cmd_string, char* params_string, char* savefile_string, char* utfile_string); +int enstrophy( int K1, int K2, int N1, int N2, double final_time, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_norm, _Complex double* u0, _Complex double* g, bool irreversible, bool keep_en_cst, double target_en, unsigned int algorithm, double print_freq, double starting_time, unsigned int nthreads, FILE* savefile, FILE* utfile, char* cmd_string, char* params_string, char* savefile_string, char* utfile_string); // 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, double final_time, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_norm, double starting_time, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int algorithm, unsigned int nthreads, FILE* savefile); +int quiet( int K1, int K2, int N1, int N2, double final_time, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_norm, double starting_time, _Complex double* u0, _Complex double* g, bool irreversible, bool keep_en_cst, double target_en, unsigned int algorithm, unsigned int nthreads, FILE* savefile); // initialize vectors for computation @@ -52,15 +52,15 @@ int copy_u( _Complex double* u, _Complex double* u0, int K1, int K2); // next time step for Irreversible/reversible Navier-Stokes equation // RK4 -int ns_step_rk4( _Complex double* u, int K1, int K2, int N1, int N2, double nu, double delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2,fft_vect ifft, _Complex double* tmp1, _Complex double *tmp2, _Complex double *tmp3, bool irreversible); +int ns_step_rk4( _Complex double* u, int K1, int K2, int N1, int N2, double nu, double delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2,fft_vect ifft, _Complex double* tmp1, _Complex double *tmp2, _Complex double *tmp3, bool irreversible, bool keep_en_cst, double target_en); // RK2 -int ns_step_rk2( _Complex double* u, int K1, int K2, int N1, int N2, double nu, double delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2,fft_vect ifft, _Complex double* tmp1, _Complex double *tmp2, bool irreversible); +int ns_step_rk2( _Complex double* u, int K1, int K2, int N1, int N2, double nu, double delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2,fft_vect ifft, _Complex double* tmp1, _Complex double *tmp2, bool irreversible, bool keep_en_cst, double target_en); // Runge-Kutta-Fehlberg -int ns_step_rkf45( _Complex double* u, double tolerance, double factor, double max_delta, unsigned int adaptive_norm, int K1, int K2, int N1, int N2, double nu, double* delta, double* next_delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double* k1, _Complex double* k2, _Complex double* k3, _Complex double* k4, _Complex double* k5, _Complex double* k6, _Complex double* tmp, bool irreversible, bool compute_k1); +int ns_step_rkf45( _Complex double* u, double tolerance, double factor, double max_delta, unsigned int adaptive_norm, int K1, int K2, int N1, int N2, double nu, double* delta, double* next_delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double* k1, _Complex double* k2, _Complex double* k3, _Complex double* k4, _Complex double* k5, _Complex double* k6, _Complex double* tmp, bool irreversible, bool keep_en_cst, double target_en, bool compute_k1); // Runge-Kutta-Dromand-Prince -int ns_step_rkdp54( _Complex double* u, double tolerance, double factor, double max_delta, unsigned int adaptive_norm, int K1, int K2, int N1, int N2, double nu, double* delta, double* next_delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double** k1, _Complex double** k2, _Complex double* k3, _Complex double* k4, _Complex double* k5, _Complex double* k6, _Complex double* tmp, bool irreversible, bool compute_k1); +int ns_step_rkdp54( _Complex double* u, double tolerance, double factor, double max_delta, unsigned int adaptive_norm, int K1, int K2, int N1, int N2, double nu, double* delta, double* next_delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double** k1, _Complex double** k2, _Complex double* k3, _Complex double* k4, _Complex double* k5, _Complex double* k6, _Complex double* tmp, bool irreversible, bool keep_en_cst, double target_en, bool compute_k1); // Runge-Kutta-Bogacki-Shampine -int ns_step_rkbs32( _Complex double* u, double tolerance, double factor, double max_delta, unsigned int adaptive_norm, int K1, int K2, int N1, int N2, double nu, double* delta, double* next_delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double** k1, _Complex double* k2, _Complex double* k3, _Complex double** k4, _Complex double* tmp, bool irreversible, bool compute_k1); +int ns_step_rkbs32( _Complex double* u, double tolerance, double factor, double max_delta, unsigned int adaptive_norm, int K1, int K2, int N1, int N2, double nu, double* delta, double* next_delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double** k1, _Complex double* k2, _Complex double* k3, _Complex double** k4, _Complex double* tmp, bool irreversible, bool keep_en_cst, double target_en, bool compute_k1); // right side of Irreversible/reversible Navier-Stokes equation int ns_rhs( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2, double nu, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft, bool irreversible);