From 814413b6e10ba90f47d7800d479ffeafcf49e67e Mon Sep 17 00:00:00 2001 From: Ian Jauslin Date: Wed, 7 Jun 2023 21:02:27 -0400 Subject: [PATCH] 1/k^3 scaling in adaptive RK --- src/navier-stokes.c | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/navier-stokes.c b/src/navier-stokes.c index 7094a3f..3471ab9 100644 --- a/src/navier-stokes.c +++ b/src/navier-stokes.c @@ -719,10 +719,11 @@ int ns_step_rkf45( for(kx=0;kx<=K1;kx++){ for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){ // 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)])); + // use the norm |u_i|/k^3 + 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)]))/pow(kx*kx+ky*ky,1.5); // 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)]); - relative+=cabs(tmp[klookup_sym(kx,ky,K2)]); + relative+=cabs(tmp[klookup_sym(kx,ky,K2)])/pow(kx*kx+ky*ky,1.5); } } @@ -816,8 +817,8 @@ int ns_step_rkbs32( for(kx=0;kx<=K1;kx++){ for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){ // difference between 5th order and 4th order - err+=cabs((*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)])); - relative+=cabs(tmp[klookup_sym(kx,ky,K2)]-u[klookup_sym(kx,ky,K2)]); + err+=cabs((*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)]))/pow(kx*kx+ky*ky,0.75); + relative+=cabs(tmp[klookup_sym(kx,ky,K2)]-u[klookup_sym(kx,ky,K2)])/pow(kx*kx+ky*ky,0.75); } } @@ -942,8 +943,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 - 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)])); - relative+=cabs(tmp[klookup_sym(kx,ky,K2)]-u[klookup_sym(kx,ky,K2)]); + // 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); } }