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

Error in enstrophy norm for rkf45 and rkbs32

Browse Source
This commit is contained in:
Ian Jauslin 2023-09-19 11:53:34 -04:00
parent d986f86243
commit f68c179038
1 changed files with 2 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/navier-stokes.c
View File

@ -815,7 +815,7 @@ int ns_step_rkf45(
err+=(kx*kx+ky*ky)*cabs2((*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)*(25./216*k1[klookup_sym(kx,ky,K2)]+1408./2565*k3[klookup_sym(kx,ky,K2)]+2197./4104*k4[klookup_sym(kx,ky,K2)]-1./5*k5[klookup_sym(kx,ky,K2)]);
sumU+=(kx*kx+ky*ky)*cabs2((*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)]));
sumU+=(kx*kx+ky*ky)*cabs2(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)]));
sumu+=(kx*kx+ky*ky)*cabs2(u[klookup_sym(kx,ky,K2)]+tmp[klookup_sym(kx,ky,K2)]);
}
}
@ -961,7 +961,7 @@ int ns_step_rkbs32(
for(kx=0;kx<=K1;kx++){
for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
err+=(kx*kx+ky*ky)*cabs2((*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)]));
sumU+=(kx*kx+ky*ky)*cabs2((*delta)*(7./24*(*k1)[klookup_sym(kx,ky,K2)]+1./4*k2[klookup_sym(kx,ky,K2)]+1./3*k3[klookup_sym(kx,ky,K2)]+1./8*(*k4)[klookup_sym(kx,ky,K2)]));
sumU+=(kx*kx+ky*ky)*cabs2(u[klookup_sym(kx,ky,K2)]+(*delta)*(7./24*(*k1)[klookup_sym(kx,ky,K2)]+1./4*k2[klookup_sym(kx,ky,K2)]+1./3*k3[klookup_sym(kx,ky,K2)]+1./8*(*k4)[klookup_sym(kx,ky,K2)]));
sumu+=(kx*kx+ky*ky)*cabs2(tmp[klookup_sym(kx,ky,K2)]);
}
}