From 051cd84a29aa0faafee6abf78c1a8a6a596882ad Mon Sep 17 00:00:00 2001 From: Ian Jauslin Date: Wed, 17 May 2023 17:41:00 -0400 Subject: [PATCH] RKDP54 --- README.md | 8 ++- src/constants.cpp | 3 +- src/main.c | 6 ++ src/navier-stokes.c | 139 +++++++++++++++++++++++++++++++++++++++++++- src/navier-stokes.h | 4 +- 5 files changed, 153 insertions(+), 7 deletions(-) diff --git a/README.md b/README.md index 97150f0..b605b22 100644 --- a/README.md +++ b/README.md @@ -93,9 +93,11 @@ should be a `;` sperated list of `key=value` pairs. The possible keys are * `random_seed` (int, default ): seed for random initialization. -* `algorithm`: `RK4` for Runge-Kutta 4, `RK2` for Runge-Kutta 2, `RKF45` for - the Runge-Kutta-Fehlberg adaptive step method, `RKBS23` for the - Runge-Kutta-Bogacki-Shampine method. +* `algorithm`: fixed step methods: `RK4` for Runge-Kutta 4, `RK2` for + Runge-Kutta 2. + adaptive step methods: `RKF45` for Runge-Kutta-Fehlberg (order + 4), `RKDP54` for Runge-Kutta-Dormand-Prince (order 5), `RKBS32` for + Runge-Kutta-Bogacki-Shampine (order 3). * `adaptive_tolerance` (double, default 1e-11): when using an adaptive step method, this is the maximal allowed relative error. diff --git a/src/constants.cpp b/src/constants.cpp index 18ed222..da593b8 100644 --- a/src/constants.cpp +++ b/src/constants.cpp @@ -34,5 +34,6 @@ limitations under the License. #define ALGORITHM_ADAPTIVE_THRESHOLD 1000 // adaptive algorithms: index > ALGORITHM_ADAPTIVE_THRESHOLD #define ALGORITHM_RKF45 1001 -#define ALGORITHM_RKBS32 1002 +#define ALGORITHM_RKDP54 1002 +#define ALGORITHM_RKBS32 1003 diff --git a/src/main.c b/src/main.c index 43b5218..d08543f 100644 --- a/src/main.c +++ b/src/main.c @@ -265,6 +265,9 @@ int print_params( case ALGORITHM_RKF45: fprintf(file,", algorithm=RKF45, tolerance=%.15e, factor=%.15e",parameters.adaptive_tolerance, parameters.adaptive_factor); break; + case ALGORITHM_RKDP54: + fprintf(file,", algorithm=RKDP54, tolerance=%.15e, factor=%.15e",parameters.adaptive_tolerance, parameters.adaptive_factor); + break; case ALGORITHM_RKBS32: fprintf(file,", algorithm=RKBS32, tolerance=%.15e, factor=%.15e",parameters.adaptive_tolerance, parameters.adaptive_factor); break; @@ -674,6 +677,9 @@ int set_parameter( else if (strcmp(rhs,"RKF45")==0){ parameters->algorithm=ALGORITHM_RKF45; } + else if (strcmp(rhs,"RKDP54")==0){ + parameters->algorithm=ALGORITHM_RKDP54; + } else if (strcmp(rhs,"RKBS32")==0){ parameters->algorithm=ALGORITHM_RKBS32; } diff --git a/src/navier-stokes.c b/src/navier-stokes.c index 12338bd..75fe2ab 100644 --- a/src/navier-stokes.c +++ b/src/navier-stokes.c @@ -90,6 +90,9 @@ int uk( ns_step_rk4(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible); } else if (algorithm==ALGORITHM_RKF45) { ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, true); + } else if (algorithm==ALGORITHM_RKDP54) { + // only compute k1 on the first step + ns_step_rkdp54(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, time==starting_time); } else if (algorithm==ALGORITHM_RKBS32) { // only compute k1 on the first step ns_step_rkbs32(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, tmp2, tmp3, &tmp4, tmp5, irreversible, time==starting_time); @@ -196,6 +199,9 @@ int enstrophy( ns_step_rk4(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible); } else if (algorithm==ALGORITHM_RKF45) { ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, true); + } else if (algorithm==ALGORITHM_RKDP54) { + // only compute k1 on the first step + ns_step_rkdp54(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, time==starting_time); } else if (algorithm==ALGORITHM_RKBS32) { // only compute k1 on the first step ns_step_rkbs32(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, tmp2, tmp3, &tmp4, tmp5, irreversible, time==starting_time); @@ -351,6 +357,9 @@ int quiet( ns_step_rk4(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible); } else if (algorithm==ALGORITHM_RKF45) { ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, true); + } else if (algorithm==ALGORITHM_RKDP54) { + // only compute k1 on the first step + ns_step_rkdp54(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, time==starting_time); } else if (algorithm==ALGORITHM_RKBS32) { // only compute k1 on the first step ns_step_rkbs32(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, tmp2, tmp3, &tmp4, tmp5, irreversible, time==starting_time); @@ -401,7 +410,7 @@ int ns_init_tmps( *tmp1=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2); *tmp2=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2); *tmp3=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2); - } else if (algorithm==ALGORITHM_RKF45){ + } else if (algorithm==ALGORITHM_RKF45 || algorithm==ALGORITHM_RKDP54){ *tmp1=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2); *tmp2=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2); *tmp3=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2); @@ -466,7 +475,7 @@ int ns_free_tmps( free(tmp1); free(tmp2); free(tmp3); - } else if (algorithm==ALGORITHM_RKF45){ + } else if (algorithm==ALGORITHM_RKF45 || algorithm==ALGORITHM_RKDP54){ free(tmp1); free(tmp2); free(tmp3); @@ -838,6 +847,132 @@ int ns_step_rkbs32( return 0; } +// next time step +// adaptive RK algorithm (Runge-Kutta-Dormand-Prince method) +int ns_step_rkdp54( + _Complex double* u, + double tolerance, + double factor, + 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, + // the pointers k1 and k2 will be exchanged at the end of the routine + _Complex double** k1, + _Complex double** k2, + _Complex double* k3, + _Complex double* k4, + _Complex double* k5, + _Complex double* k6, + _Complex double* tmp, + bool irreversible, + // whether to compute k1 + bool compute_k1 +){ + int kx,ky; + double err,relative; + + // k1: u(t) + // only compute it if it is the first step (otherwise, it has already been computed due to the First Same As Last property) + if(compute_k1){ + ns_rhs(*k1, u, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible); + } + + // k2 : u(t+1/5*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)/5*(*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/10*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./40*(*k1)[klookup_sym(kx,ky,K2)]+9./40*(*k2)[klookup_sym(kx,ky,K2)]); + } + } + ns_rhs(k3, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible); + + // k4 : u(t+4/5*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)*(44./45*(*k1)[klookup_sym(kx,ky,K2)]-56./15*(*k2)[klookup_sym(kx,ky,K2)]+32./9*k3[klookup_sym(kx,ky,K2)]); + } + } + ns_rhs(k4, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible); + + // k5 : u(t+8/9*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)*(19372./6561*(*k1)[klookup_sym(kx,ky,K2)]-25360./2187*(*k2)[klookup_sym(kx,ky,K2)]+64448./6561*k3[klookup_sym(kx,ky,K2)]-212./729*k4[klookup_sym(kx,ky,K2)]); + } + } + ns_rhs(k5, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible); + + // k6 : u(t+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)*(9017./3168*(*k1)[klookup_sym(kx,ky,K2)]-355./33*(*k2)[klookup_sym(kx,ky,K2)]+46732./5247*k3[klookup_sym(kx,ky,K2)]+49./176*k4[klookup_sym(kx,ky,K2)]-5103./18656*k5[klookup_sym(kx,ky,K2)]); + } + } + ns_rhs(k6, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible); + + // k7 : u(t+delta) + for(kx=0;kx<=K1;kx++){ + for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){ + // tmp computed here is the step + tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(35./384*(*k1)[klookup_sym(kx,ky,K2)]+500./1113*k3[klookup_sym(kx,ky,K2)]+125./192*k4[klookup_sym(kx,ky,K2)]-2187./6784*k5[klookup_sym(kx,ky,K2)]+11./84*k6[klookup_sym(kx,ky,K2)]); + } + } + // store in k2, which is not needed anymore + ns_rhs(*k2, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible); + + // compute error + err=0; + relative=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)*(-71./57600*(*k1)[klookup_sym(kx,ky,K2)]+71./16695*k3[klookup_sym(kx,ky,K2)]-71./1920*k4[klookup_sym(kx,ky,K2)]+17253./339200*k5[klookup_sym(kx,ky,K2)]-22./525*k6[klookup_sym(kx,ky,K2)]+1./40*(*k2)[klookup_sym(kx,ky,K2)])); + relative+=cabs(tmp[klookup_sym(kx,ky,K2)]-u[klookup_sym(kx,ky,K2)]); + } + } + + // compare relative error with tolerance + if(err0 ? -K2 : 1);ky<=K2;ky++){ + u[klookup_sym(kx,ky,K2)]=tmp[klookup_sym(kx,ky,K2)]; + } + } + // next delta to use in future steps + *next_delta=(*delta)*pow(relative*tolerance/err,1./5); + + // k1 in the next step will be this k4 (first same as last) + tmp=*k1; + *k1=*k2; + *k2=tmp; + } + // 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,K1,K2,N1,N2,nu,delta,next_delta,L,g,fft1,fft2,ifft,k1,k2,k3,k4,k5,k6,tmp,irreversible,false); + } + + return 0; +} + // 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 ec35052..8a57a53 100644 --- a/src/navier-stokes.h +++ b/src/navier-stokes.h @@ -55,8 +55,10 @@ 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) +// Runge-Kutta-Fehlberg int ns_step_rkf45( _Complex double* u, double tolerance, double factor, 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); +// Runge-Kutta-Dromand-Prince +int ns_step_rkdp54( _Complex double* u, double tolerance, double factor, 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); // Runge-Kutta-Bogacki-Shampine int ns_step_rkbs32( _Complex double* u, double tolerance, double factor, 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);