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

Use the enstrophy-norm in rkdp54

Browse Source
This commit is contained in:
Ian Jauslin 2023-06-13 17:30:51 -04:00
parent 814413b6e1
commit 4ffbe1e978
1 changed files with 5 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
src/navier-stokes.c
View File

@ -22,6 +22,8 @@ limitations under the License.
#include <stdlib.h>
#include <string.h>
#define CABS2(x) ((__real__ x) * (__real__ x) + (__imag__ x) * (__imag__ x))
// compute solution as a function of time
int uk(
int K1,
@ -943,9 +945,9 @@ int ns_step_rkdp54(
for(kx=0;kx<=K1;kx++){
for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
// difference between 5th order and 4th order
// use the norm |u_i|/k^3
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)]))/pow(kx*kx+ky*ky,1.5);
relative+=cabs(tmp[klookup_sym(kx,ky,K2)]-u[klookup_sym(kx,ky,K2)])/pow(kx*kx+ky*ky,1.5);
// use the norm |u_k|^2k^2 (to get a bound on the error of the enstrophy)
err+=(kx*kx+ky*ky)*CABS2((*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+=(kx*kx+ky*ky)*(CABS2(tmp[klookup_sym(kx,ky,K2)]-u[klookup_sym(kx,ky,K2)]));
}
}