software/Nstrophy
1
Fork 0
Code Issues Pull Requests Releases Activity

RKBS23

Browse Source
This commit is contained in:
Ian Jauslin 2023-05-17 16:47:15 -04:00
parent c9312c6d16
commit d069e8c94e
5 changed files with 131 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
README.md
View File

@ -93,8 +93,9 @@ 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, `RK2` for Runge-Kutta 2, or
`RKF45` for the Runge-Kutta-Fehlberg adaptive step method.
* `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.
* `adaptive_tolerance` (double, default 1e-11): when using an adaptive step
method, this is the maximal allowed relative error.

1
src/constants.cpp
View File

@ -34,4 +34,5 @@ limitations under the License.
#define ALGORITHM_ADAPTIVE_THRESHOLD 1000
// adaptive algorithms: index > ALGORITHM_ADAPTIVE_THRESHOLD
#define ALGORITHM_RKF45 1001
#define ALGORITHM_RKBS23 1002

6
src/main.c
View File

@ -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_RKBS23:
fprintf(file,", algorithm=RKBS23, tolerance=%.15e, factor=%.15e",parameters.adaptive_tolerance, parameters.adaptive_factor);
break;
default:
fprintf(file,", algorithm=RK4");
break;
@ -671,6 +674,9 @@ int set_parameter(
else if (strcmp(rhs,"RKF45")==0){
parameters->algorithm=ALGORITHM_RKF45;
}
else if (strcmp(rhs,"RKBS23")==0){
parameters->algorithm=ALGORITHM_RKBS23;
}
else{
fprintf(stderr, "error: unrecognized algorithm '%s'\n",rhs);
return(-1);

120
src/navier-stokes.c
View File

@ -90,6 +90,8 @@ 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);
} else if (algorithm==ALGORITHM_RKBS23) {
ns_step_rkbs23(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);
} else {
fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
}
@ -193,6 +195,8 @@ 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);
} else if (algorithm==ALGORITHM_RKBS23) {
ns_step_rkbs23(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);
} else {
fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
}
@ -216,7 +220,7 @@ int enstrophy(
avg_en_x_a*=print_freq/(time-prevtime);
// print to stderr so user can follow along
if(algorithm==ALGORITHM_RKF45){
if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){
fprintf(stderr,"% .8e % .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",time, avg_a, avg_en, avg_en_x_a, alpha, enstrophy, alpha*enstrophy, step);
printf("% .15e % .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",time, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy, step);
} else {
@ -345,6 +349,8 @@ 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);
} else if (algorithm==ALGORITHM_RKBS23) {
ns_step_rkbs23(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);
} else {
fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
}
@ -400,6 +406,12 @@ int ns_init_tmps(
*tmp5=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
*tmp6=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
*tmp7=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
} else if (algorithm==ALGORITHM_RKBS23){
*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);
*tmp4=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
*tmp5=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
} else {
fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
};
@ -459,6 +471,12 @@ int ns_free_tmps(
free(tmp5);
free(tmp6);
free(tmp7);
} else if (algorithm==ALGORITHM_RKBS23){
free(tmp1);
free(tmp2);
free(tmp3);
free(tmp4);
free(tmp5);
} else {
fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
};
@ -712,6 +730,106 @@ int ns_step_rkf45(
return 0;
}
// next time step
// adaptive RK algorithm (Runge-Kutta-Bogacki-Shampine method)
int ns_step_rkbs23(
_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 k4 will be exchanged at the end of the routine
_Complex double** k1,
_Complex double* k2,
_Complex double* k3,
_Complex double** k4,
_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/4*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)/2*(*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/4*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./4*k2[klookup_sym(kx,ky,K2)]);
}
}
ns_rhs(k3, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
// k4 : u(t+delta)
// tmp cpmputed here is the next step
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)*(2./9*(*k1)[klookup_sym(kx,ky,K2)]+1./3*k2[klookup_sym(kx,ky,K2)]+4./9*k3[klookup_sym(kx,ky,K2)]);
}
}
ns_rhs(*k4, 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)*(5./72*(*k1)[klookup_sym(kx,ky,K2)]-1./12*k2[klookup_sym(kx,ky,K2)]-1./9*k3[klookup_sym(kx,ky,K2)]+1./8*(*k4)[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(err<relative*tolerance){
// add to output
for(kx=0;kx<=K1;kx++){
for(ky=(kx>0 ? -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./3);
// k1 in the next step will be this k4 (first same as last)
tmp=*k1;
*k1=*k4;
*k4=tmp;
}
// 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_rkbs23(u,tolerance,factor,K1,K2,N1,N2,nu,delta,next_delta,L,g,fft1,fft2,ifft,k1,k2,k3,k4,tmp,irreversible,false);
}
return 0;
}
// right side of Irreversible/Reversible Navier-Stokes equation
int ns_rhs(
_Complex double* out,

2
src/navier-stokes.h
View File

@ -57,6 +57,8 @@ int ns_step_rk4( _Complex double* u, int K1, int K2, int N1, int N2, double nu,
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)
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);
// Runge-Kutta-Bogacki-Shampine
int ns_step_rkbs23( _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);
// 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);