From 5db3680b31e25eda7517aff37a1c3ffa7fa7cc86 Mon Sep 17 00:00:00 2001 From: Ian Jauslin Date: Mon, 15 May 2023 20:29:06 -0400 Subject: [PATCH] RKF45 algorithm --- README.md | 13 ++- src/constants.cpp | 5 +- src/main.c | 86 ++++++++++---- src/navier-stokes.c | 271 ++++++++++++++++++++++++++++++++++++++------ src/navier-stokes.h | 12 +- 5 files changed, 325 insertions(+), 62 deletions(-) diff --git a/README.md b/README.md index f7c01fe..fa2f6f0 100644 --- a/README.md +++ b/README.md @@ -68,7 +68,9 @@ should be a `;` sperated list of `key=value` pairs. The possible keys are * `print_freq` (long int, default 1000): only print every `print_freq` steps. -* `starting_time` (long int, default 0): start the computation at this step. +* `starting_step` (long int, default 0): start the computation at this step. + +* `starting_time` (double, default starting_time*delta): start the computation at this time. * `driving`: either `zero` for no driving, `test` (default) for a testing driving force or `file:` or `file_txt:` to read the @@ -91,7 +93,14 @@ should be a `;` sperated list of `key=value` pairs. The possible keys are * `random_seed` (int, default ): seed for random initialization. -* `algorithm`: either `RK4` for Runge-Kutta 4, or `RK2` for Runge-Kutta 2. +* `algorithm`: either `RK4` for Runge-Kutta 4, `RK2` for Runge-Kutta 2, or + `RKF45` for the Runge-Kutta-Fehlberg adaptive step method. + +* `adaptive_tolerance` (double, default 1e-11): when using an adaptive step + method, this is the maximal allowed relative error. + +* `adaptive_factor` (double, default 0.9): when using the RKF45 method, the + step gets adjusted by `factor*delta*(tolerance/error)^(1/5)`. # Interrupting/resuming the computation diff --git a/src/constants.cpp b/src/constants.cpp index bc49a04..8bd85ed 100644 --- a/src/constants.cpp +++ b/src/constants.cpp @@ -29,6 +29,9 @@ limitations under the License. #define INIT_FILE 3 #define INIT_FILE_TXT 4 -#define ALGORITHM_RK4 4 +#define ALGORITHM_RK4 1 #define ALGORITHM_RK2 2 +#define ALGORITHM_ADAPTIVE_THRESHOLD 1000 +// adaptive algorithms: index > ALGORITHM_ADAPTIVE_THRESHOLD +#define ALGORITHM_RKF45 1001 diff --git a/src/main.c b/src/main.c index e32569c..bc4cc0c 100644 --- a/src/main.c +++ b/src/main.c @@ -46,9 +46,12 @@ typedef struct nstrophy_parameters { double nu; double delta; double L; + double adaptive_tolerance; + double adaptive_factor; uint64_t print_freq; int seed; - uint64_t starting_time; + uint64_t starting_step; + double starting_time; unsigned int driving; unsigned int init; unsigned int algorithm; @@ -64,12 +67,12 @@ int print_params(nstrophy_parameters parameters, char* initfile_str, char* drivi // read command line arguments int read_args(int argc, const char* argv[], char** params, unsigned int* command, unsigned int* nthreads, char** savefile_str); int read_params(char* param_str, nstrophy_parameters* parameters, char** initfile_str, char** drivingfile_str); -int set_parameter(char* lhs, char* rhs, nstrophy_parameters* parameters, bool* setN1, bool* setN2, char** initfile_str, char** drivingfile_str); +int set_parameter(char* lhs, char* rhs, nstrophy_parameters* parameters, bool* setN1, bool* setN2, bool* set_starting_time, char** initfile_str, char** drivingfile_str); // set driving force _Complex double* set_driving(nstrophy_parameters parameters); // set initial condition -_Complex double* set_init(nstrophy_parameters parameters); +_Complex double* set_init(nstrophy_parameters* parameters); // signal handler void sig_handler (int signo); @@ -134,7 +137,7 @@ int main ( // set driving force g=set_driving(parameters); // set initial condition - u0=set_init(parameters); + u0=set_init(¶meters); // close initfile (do this early, so that it is possible to use the same file for init and save) if (parameters.initfile!=NULL){ @@ -169,15 +172,15 @@ int main ( // run command if (command==COMMAND_UK){ - uk(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.irreversible, parameters.algorithm, parameters.print_freq, parameters.starting_time, nthreads, savefile); + uk(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, u0, g, parameters.irreversible, parameters.algorithm, parameters.print_freq, parameters.starting_step, parameters.starting_time, nthreads, savefile); } else if(command==COMMAND_ENSTROPHY){ // register signal handler to handle aborts signal(SIGINT, sig_handler); - enstrophy(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.irreversible, parameters.algorithm, parameters.print_freq, parameters.starting_time, nthreads, savefile, (char*)argv[0], param_str, savefile_str); + enstrophy(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, u0, g, parameters.irreversible, parameters.algorithm, parameters.print_freq, parameters.starting_step, parameters.starting_time, nthreads, savefile, (char*)argv[0], param_str, savefile_str); } else if(command==COMMAND_QUIET){ - quiet(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, parameters.starting_time, u0, g, parameters.irreversible, parameters.algorithm, nthreads, savefile); + quiet(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.starting_step, u0, g, parameters.irreversible, parameters.algorithm, nthreads, savefile); } else if(command==0){ fprintf(stderr, "error: no command specified\n"); @@ -261,6 +264,9 @@ int print_params( case ALGORITHM_RK2: fprintf(file,", algorithm=RK2"); break; + case ALGORITHM_RKF45: + fprintf(file,", algorithm=RKF45, tolerance=%.15e, factor=%.15e",parameters.adaptive_tolerance, parameters.adaptive_factor); + break; default: fprintf(file,", algorithm=RK4"); break; @@ -371,6 +377,7 @@ int read_params( // whether N was set explicitly bool setN1=false; bool setN2=false; + bool set_starting_time=false; // whether lhs (false is rhs) bool lhs=true; @@ -384,9 +391,12 @@ int read_params( //nu=2^-11 parameters->nu=0.00048828125; parameters->L=2*M_PI; + parameters->adaptive_tolerance=1e-11; + parameters->adaptive_factor=0.9; parameters->nsteps=10000000; parameters->print_freq=1000; parameters->starting_time=0; + parameters->starting_step=0; parameters->seed=17; parameters->driving=DRIVING_TEST; parameters->drivingfile=NULL; @@ -413,7 +423,7 @@ int read_params( break; case ';': //set parameter - ret=set_parameter(buffer_lhs, buffer_rhs, parameters, &setN1, &setN2, initfile_str, drivingfile_str); + ret=set_parameter(buffer_lhs, buffer_rhs, parameters, &setN1, &setN2, &set_starting_time, initfile_str, drivingfile_str); if(ret<0){ return ret; } @@ -440,7 +450,7 @@ int read_params( // set last param if (*param_str!='\0'){ - ret=set_parameter(buffer_lhs, buffer_rhs, parameters, &setN1, &setN2, initfile_str, drivingfile_str); + ret=set_parameter(buffer_lhs, buffer_rhs, parameters, &setN1, &setN2, &set_starting_time, initfile_str, drivingfile_str); if(ret<0){ return ret; } @@ -459,6 +469,10 @@ int read_params( parameters->N2=smallest_pow2(3*(parameters->K2)); } + // if starting_time not set explicitly and algorithm is not adaptive, set it to delta*starting_step + if (!set_starting_time && parameters->algorithmstarting_time=parameters->delta*parameters->starting_step; + } return(0); } @@ -470,6 +484,7 @@ int set_parameter( nstrophy_parameters* parameters, bool* setN1, bool* setN2, + bool* set_starting_time, char** initfile_str, char** drivingfile_str ){ @@ -570,6 +585,20 @@ int set_parameter( return(-1); } } + else if (strcmp(lhs,"adaptive_tolerance")==0){ + ret=sscanf(rhs,"%lf",&(parameters->adaptive_tolerance)); + if(ret!=1){ + fprintf(stderr, "error: parameter 'adaptive_tolerance' should be a double\n got '%s'\n",rhs); + return(-1); + } + } + else if (strcmp(lhs,"adaptive_factor")==0){ + ret=sscanf(rhs,"%lf",&(parameters->adaptive_factor)); + if(ret!=1){ + fprintf(stderr, "error: parameter 'adaptive_factor' should be a double\n got '%s'\n",rhs); + return(-1); + } + } else if (strcmp(lhs,"print_freq")==0){ ret=sscanf(rhs,"%lu",&(parameters->print_freq)); if(ret!=1){ @@ -584,13 +613,21 @@ int set_parameter( return(-1); } } - else if (strcmp(lhs,"starting_time")==0){ - ret=sscanf(rhs,"%lu",&(parameters->starting_time)); + else if (strcmp(lhs,"starting_step")==0){ + ret=sscanf(rhs,"%lu",&(parameters->starting_step)); if(ret!=1){ - fprintf(stderr, "error: parameter 'starting_time' should be an unsigned integer\n got '%s'\n",rhs); + fprintf(stderr, "error: parameter 'starting_step' should be a long unsigned integer\n got '%s'\n",rhs); return(-1); } } + else if (strcmp(lhs,"starting_time")==0){ + ret=sscanf(rhs,"%lf",&(parameters->starting_time)); + if(ret!=1){ + fprintf(stderr, "error: parameter 'starting_time' should be a double\n got '%s'\n",rhs); + return(-1); + } + *set_starting_time=true; + } else if (strcmp(lhs,"driving")==0){ if (strcmp(rhs,"zero")==0){ parameters->driving=DRIVING_ZERO; @@ -648,6 +685,9 @@ int set_parameter( else if (strcmp(rhs,"RK2")==0){ parameters->algorithm=ALGORITHM_RK2; } + else if (strcmp(rhs,"RKF45")==0){ + parameters->algorithm=ALGORITHM_RKF45; + } else{ fprintf(stderr, "error: unrecognized algorithm '%s'\n",rhs); return(-1); @@ -692,29 +732,35 @@ _Complex double* set_driving( // set initial condition _Complex double* set_init( - nstrophy_parameters parameters + nstrophy_parameters* parameters ){ - _Complex double* u0=calloc(sizeof(_Complex double),parameters.K1*(2*parameters.K2+1)+parameters.K2); + _Complex double* u0=calloc(sizeof(_Complex double),parameters->K1*(2*parameters->K2+1)+parameters->K2); - switch(parameters.init){ + switch(parameters->init){ case INIT_RANDOM: - init_random(u0, parameters.init_en, parameters.K1, parameters.K2, parameters.L, parameters.seed, parameters.irreversible); + init_random(u0, parameters->init_en, parameters->K1, parameters->K2, parameters->L, parameters->seed, parameters->irreversible); break; case INIT_GAUSSIAN: - init_gaussian(u0, parameters.init_en, parameters.K1, parameters.K2, parameters.L, parameters.irreversible); + init_gaussian(u0, parameters->init_en, parameters->K1, parameters->K2, parameters->L, parameters->irreversible); break; case INIT_FILE: - init_file_bin(u0, parameters.K1, parameters.K2, parameters.initfile); + init_file_bin(u0, parameters->K1, parameters->K2, parameters->initfile); + if(parameters->algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){ + // read delta + fread(&(parameters->delta), sizeof(double), 1, parameters->initfile); + // read start time + fread(&(parameters->starting_time), sizeof(double), 1, parameters->initfile); + } break; case INIT_FILE_TXT: - init_file_txt(u0, parameters.K1, parameters.K2, parameters.initfile); + init_file_txt(u0, parameters->K1, parameters->K2, parameters->initfile); break; default: - init_gaussian(u0, parameters.init_en, parameters.K1, parameters.K2, parameters.L, parameters.irreversible); + init_gaussian(u0, parameters->init_en, parameters->K1, parameters->K2, parameters->L, parameters->irreversible); break; } diff --git a/src/navier-stokes.c b/src/navier-stokes.c index b5fe7c7..8594e5c 100644 --- a/src/navier-stokes.c +++ b/src/navier-stokes.c @@ -32,26 +32,34 @@ int uk( double nu, double delta, double L, + double adaptive_tolerance, + double adaptive_factor, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int algorithm, uint64_t print_freq, - uint64_t starting_time, + uint64_t starting_step, + double starting_time, unsigned int nthreads, FILE* savefile ){ + double time=starting_time; _Complex double* u; _Complex double* tmp1; _Complex double* tmp2; _Complex double* tmp3; + _Complex double* tmp4; + _Complex double* tmp5; + _Complex double* tmp6; + _Complex double* tmp7; uint64_t t; fft_vect fft1; fft_vect fft2; fft_vect ifft; int kx,ky; - ns_init_tmps(&u, &tmp1, &tmp2, &tmp3, &fft1, &fft2, &ifft, K1, K2, N1, N2, nthreads); + ns_init_tmps(&u, &tmp1, &tmp2, &tmp3, &tmp4, &tmp5, &tmp6, &tmp7, &fft1, &fft2, &ifft, K1, K2, N1, N2, nthreads, algorithm); // copy initial condition copy_u(u, u0, K1, K2); @@ -68,16 +76,22 @@ int uk( } // iterate - for(t=starting_time;nsteps==0 || tstarting_time && t%print_freq==0){ - fprintf(stderr,"%lu % .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",t,t*delta, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy); - printf("%8lu % .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",t,t*delta, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy); + if(t>starting_step && t%print_freq==0){ + // print to stderr so user can follow along + if(algorithm==ALGORITHM_RKF45){ + fprintf(stderr,"%lu % .8e % .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",t,time, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy, delta); + } else { + fprintf(stderr,"%lu % .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",t,time, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy); + } + // print to stdout + printf("%8lu % .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",t,time, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy); } // catch abort signal @@ -189,11 +223,24 @@ int enstrophy( if(params_string!=NULL) { char* params=calloc(sizeof(char), strlen(params_string)+1); strcpy(params, params_string); + remove_entry(params, "starting_step"); remove_entry(params, "starting_time"); remove_entry(params, "init"); remove_entry(params, "nsteps"); - fprintf(savefile," -p \"%s;starting_time=%lu;nsteps=%lu;init=file:%s\"", params, t+1, (nsteps+starting_time < t+1 ? 0 : nsteps+starting_time-t-1), savefile_string); + if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){ + remove_entry(params, "delta"); + } + fprintf(savefile," -p \"%s;starting_step=%lu;nsteps=%lu;init=file:%s", params, t+1, (nsteps+starting_step < t+1 ? 0 : nsteps+starting_step-t-1), savefile_string); free(params); + // write delta if adaptive, and not writing binary + if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD && (savefile==stderr || savefile==stdout)){ + fprintf(savefile,";delta=%.15e;starting_time=%.15e", delta, starting_time); + } + // write starting_time if not adaptive + if(algorithmALGORITHM_ADAPTIVE_THRESHOLD){ + // first binary entry: delta + fwrite(&delta, sizeof(double), 1, savefile); + // ssecond binary entry: starting time + fwrite(&starting_time, sizeof(double), 1, savefile); + } } } - ns_free_tmps(u, tmp1, tmp2, tmp3, fft1, fft2, ifft); + ns_free_tmps(u, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, fft1, fft2, ifft, algorithm); return(0); } @@ -219,7 +273,9 @@ int quiet( double nu, double delta, double L, - uint64_t starting_time, + double adaptive_tolerance, + double adaptive_factor, + uint64_t starting_step, _Complex double* u0, _Complex double* g, bool irreversible, @@ -231,38 +287,50 @@ int quiet( _Complex double* tmp1; _Complex double* tmp2; _Complex double* tmp3; + _Complex double* tmp4; + _Complex double* tmp5; + _Complex double* tmp6; + _Complex double* tmp7; uint64_t t; fft_vect fft1; fft_vect fft2; fft_vect ifft; - ns_init_tmps(&u, &tmp1, &tmp2, &tmp3, &fft1, &fft2, &ifft, K1, K2, N1, N2, nthreads); + ns_init_tmps(&u, &tmp1, &tmp2, &tmp3, &tmp4, &tmp5, &tmp6, &tmp7, &fft1, &fft2, &ifft, K1, K2, N1, N2, nthreads, algorithm); // copy initial condition copy_u(u, u0, K1, K2); // iterate - for(t=starting_time;nsteps==0 || t0 ? -K2 : 1);ky<=K2;ky++){ + tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta/4*k1[klookup_sym(kx,ky,K2)]; + } + } + ns_rhs(k2, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible); + + // k3 : u(t+3/8*delta) + for(kx=0;kx<=K1;kx++){ + for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){ + tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta*(3./32*k1[klookup_sym(kx,ky,K2)]+9./32*k2[klookup_sym(kx,ky,K2)]); + } + } + ns_rhs(k3, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible); + + // k4 : u(t+12./13*delta) + for(kx=0;kx<=K1;kx++){ + for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){ + tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta*(1932./2197*k1[klookup_sym(kx,ky,K2)]-7200./2197*k2[klookup_sym(kx,ky,K2)]+7296./2197*k3[klookup_sym(kx,ky,K2)]); + } + } + ns_rhs(k4, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible); + + // k5 : u(t+1*delta) + for(kx=0;kx<=K1;kx++){ + for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){ + tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta*(439./216*k1[klookup_sym(kx,ky,K2)]-8*k2[klookup_sym(kx,ky,K2)]+3680./513*k3[klookup_sym(kx,ky,K2)]-845./4104*k4[klookup_sym(kx,ky,K2)]); + } + } + ns_rhs(k5, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible); + + // k6 : u(t+1./2*delta) + for(kx=0;kx<=K1;kx++){ + for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){ + tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta*(-8./27*k1[klookup_sym(kx,ky,K2)]+2*k2[klookup_sym(kx,ky,K2)]-3544./2565*k3[klookup_sym(kx,ky,K2)]+1859./4104*k4[klookup_sym(kx,ky,K2)]-11./40*k5[klookup_sym(kx,ky,K2)]); + } + } + ns_rhs(k6, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible); + + // compute error + err=0; + for(kx=0;kx<=K1;kx++){ + for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){ + // difference between 5th order and 4th order + err+=cabs(delta*(1./360*k1[klookup_sym(kx,ky,K2)]-128./4275*k3[klookup_sym(kx,ky,K2)]-2197./75240*k4[klookup_sym(kx,ky,K2)]+1./50*k5[klookup_sym(kx,ky,K2)]+2./55*k6[klookup_sym(kx,ky,K2)])); + } + } + + // new delta + delta=factor*delta*pow(tolerance/err,0.2); + // compare relative error with tolerance + if(err0 ? -K2 : 1);ky<=K2;ky++){ + u[klookup_sym(kx,ky,K2)]+=delta*(16./135*k1[klookup_sym(kx,ky,K2)]+6656./12825*k3[klookup_sym(kx,ky,K2)]+28561./56430*k4[klookup_sym(kx,ky,K2)]-9./50*k5[klookup_sym(kx,ky,K2)]+2./55*k6[klookup_sym(kx,ky,K2)]); + } + } + return delta; + } + // error too big: repeat with smaller step + else{ + return(ns_step_rkf45(u,tolerance,factor,K1,K2,N1,N2,nu,delta,L,g,fft1,fft2,ifft,k1,k2,k3,k4,k5,k6,tmp,irreversible)); + } +} + // right side of Irreversible/Reversible Navier-Stokes equation int ns_rhs( _Complex double* out, diff --git a/src/navier-stokes.h b/src/navier-stokes.h index aa02d8a..331ece6 100644 --- a/src/navier-stokes.h +++ b/src/navier-stokes.h @@ -34,19 +34,19 @@ typedef struct fft_vects { } fft_vect; // compute u_k -int uk( int K1, int K2, int N1, int N2, uint64_t nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int algorithm, uint64_t print_freq, uint64_t starting_time, unsigned int nthreadsl, FILE* savefile); +int uk( int K1, int K2, int N1, int N2, uint64_t nsteps, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int algorithm, uint64_t print_freq, uint64_t starting_step, double starting_time, unsigned int nthreadsl, FILE* savefile); // compute enstrophy and alpha -int enstrophy( int K1, int K2, int N1, int N2, uint64_t nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int algorithm, uint64_t print_freq, uint64_t starting_time, unsigned int nthreads, FILE* savefile, char* cmd_string, char* params_string, char* savefile_string); +int enstrophy( int K1, int K2, int N1, int N2, uint64_t nsteps, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int algorithm, uint64_t print_freq, uint64_t starting_step, double starting_time, unsigned int nthreads, FILE* savefile, char* cmd_string, char* params_string, char* savefile_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, uint64_t nsteps, double nu, double delta, double L, uint64_t 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, uint64_t nsteps, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, uint64_t starting_step, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int algorithm, unsigned int nthreads, FILE* savefile); // initialize vectors for computation -int ns_init_tmps( _Complex double **u, _Complex double ** tmp1, _Complex double **tmp2, _Complex double **tmp3, fft_vect* fft1, fft_vect *fft2, fft_vect *ifft, int K1, int K2, int N1, int N2, unsigned int nthreads); +int ns_init_tmps( _Complex double **u, _Complex double ** tmp1, _Complex double **tmp2, _Complex double **tmp3, _Complex double **tmp4, _Complex double **tmp5, _Complex double **tmp6, _Complex double **tmp7, fft_vect* fft1, fft_vect *fft2, fft_vect *ifft, int K1, int K2, int N1, int N2, unsigned int nthreads, unsigned int algorithm); // 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); +int ns_free_tmps( _Complex double* u, _Complex double* tmp1, _Complex double *tmp2,_Complex double *tmp3, _Complex double *tmp4, _Complex double *tmp5, _Complex double *tmp6, _Complex double *tmp7, fft_vect fft1, fft_vect fft2, fft_vect ifft, unsigned int algorithm); // copy u0 to u int copy_u( _Complex double* u, _Complex double* u0, int K1, int K2); @@ -56,6 +56,8 @@ int copy_u( _Complex double* u, _Complex double* u0, int K1, int K2); 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); // 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); +// adaptive RK algorithm (Runge-Kutta-Fehlberg) +double ns_step_rkf45( _Complex double* u, double tolerance, double factor, 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* k1, _Complex double* k2, _Complex double* k3, _Complex double* k4, _Complex double* k5, _Complex double* k6, _Complex double* tmp, bool irreversible); // 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);