diff --git a/src/navier-stokes.c b/src/navier-stokes.c index 3471ab9..186864e 100644 --- a/src/navier-stokes.c +++ b/src/navier-stokes.c @@ -22,6 +22,8 @@ limitations under the License. #include #include +#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)])); } }