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

Use a relative tolerance for adaptive step size

Browse Source
This commit is contained in:
Ian Jauslin 2023-05-16 16:49:50 -04:00
parent 949b3897c4
commit cf8cdfe12c
1 changed files with 8 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
src/navier-stokes.c
View File

@ -615,7 +615,7 @@ double ns_step_rkf45(
bool irreversible bool irreversible
){ ){
int kx,ky; int kx,ky;
double err; double err,relative;
// k1: u(t) // k1: u(t)
ns_rhs(k1, u, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible); ns_rhs(k1, u, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
@ -662,21 +662,25 @@ double ns_step_rkf45(
// compute error // compute error
err=0; err=0;
relative=0;
for(kx=0;kx<=K1;kx++){ for(kx=0;kx<=K1;kx++){
for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){ for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
// difference between 5th order and 4th order // 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)])); 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)]));
// next step
tmp[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)]);
relative+=cabs(tmp[klookup_sym(kx,ky,K2)]);
} }
} }
// new delta // new delta
delta=factor*delta*pow(tolerance/err,0.2); delta=factor*delta*pow(relative*tolerance/err,0.2);
// compare relative error with tolerance // compare relative error with tolerance
if(err<tolerance){ if(err<relative*tolerance){
// add to output // add to output
for(kx=0;kx<=K1;kx++){ for(kx=0;kx<=K1;kx++){
for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){ for(ky=(kx>0 ? -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)]); u[klookup_sym(kx,ky,K2)]+=tmp[klookup_sym(kx,ky,K2)];
} }
} }
return delta; return delta;