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

RKF45: no not recompute k1 when decreasing step

Browse Source
This commit is contained in:
Ian Jauslin
2023-05-17 16:50:11 -04:00
parent d069e8c94e
commit f57cc317ac
2 changed files with 15 additions and 7 deletions
Download Patch File Download Diff File Expand all files Collapse all files

20
src/navier-stokes.c
View File

@@ -89,8 +89,9 @@ int uk(
} else if (algorithm==ALGORITHM_RK4) {
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);
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_RKBS23) {
// only compute k1 on the first step
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);
@@ -194,8 +195,9 @@ int enstrophy(
} else if (algorithm==ALGORITHM_RK4) {
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);
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_RKBS23) {
// only compute k1 on the first step
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);
@@ -348,8 +350,9 @@ int quiet(
} else if (algorithm==ALGORITHM_RK4) {
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);
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_RKBS23) {
// only compute k1 on the first step
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);
@@ -649,13 +652,17 @@ int ns_step_rkf45(
_Complex double* k5,
_Complex double* k6,
_Complex double* tmp,
bool irreversible
bool irreversible,
// whether to compute k1 or leave it as is
bool compute_k1
){
int kx,ky;
double err,relative;
// k1: u(t)
ns_rhs(k1, u, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
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++){
@@ -724,7 +731,8 @@ int ns_step_rkf45(
// error too big: repeat with smaller step
else{
*delta=factor*(*delta)*pow(relative*tolerance/err,0.2);
ns_step_rkf45(u,tolerance,factor,K1,K2,N1,N2,nu,delta,next_delta,L,g,fft1,fft2,ifft,k1,k2,k3,k4,k5,k6,tmp,irreversible);
// do not recompute k1
ns_step_rkf45(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;