From aeec3ff17b39a20f9a038f939b079692fa445b3a Mon Sep 17 00:00:00 2001 From: Ian Jauslin Date: Wed, 14 Jun 2023 00:14:17 -0400 Subject: [PATCH] max_delta and adaptive norm in rkf45 and rkbs32 --- src/navier-stokes.c | 140 ++++++++++++++++++++++++++++++++++---------- src/navier-stokes.h | 4 +- 2 files changed, 112 insertions(+), 32 deletions(-) diff --git a/src/navier-stokes.c b/src/navier-stokes.c index ce699f6..b33d3ba 100644 --- a/src/navier-stokes.c +++ b/src/navier-stokes.c @@ -92,13 +92,13 @@ int uk( } else if (algorithm==ALGORITHM_RK4) { ns_step_rk4(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible); } else if (algorithm==ALGORITHM_RKF45) { - ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, true); + ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, true); } else if (algorithm==ALGORITHM_RKDP54) { // only compute k1 on the first step ns_step_rkdp54(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, time==starting_time); } else if (algorithm==ALGORITHM_RKBS32) { // only compute k1 on the first step - ns_step_rkbs32(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, tmp2, tmp3, &tmp4, tmp5, irreversible, time==starting_time); + ns_step_rkbs32(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, tmp2, tmp3, &tmp4, tmp5, irreversible, time==starting_time); } else { fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm); } @@ -203,13 +203,13 @@ int enstrophy( } else if (algorithm==ALGORITHM_RK4) { ns_step_rk4(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible); } else if (algorithm==ALGORITHM_RKF45) { - ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, true); + ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, true); } else if (algorithm==ALGORITHM_RKDP54) { // only compute k1 on the first step ns_step_rkdp54(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, time==starting_time); } else if (algorithm==ALGORITHM_RKBS32) { // only compute k1 on the first step - ns_step_rkbs32(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, tmp2, tmp3, &tmp4, tmp5, irreversible, time==starting_time); + ns_step_rkbs32(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, tmp2, tmp3, &tmp4, tmp5, irreversible, time==starting_time); } else { fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm); } @@ -363,13 +363,13 @@ int quiet( } else if (algorithm==ALGORITHM_RK4) { ns_step_rk4(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible); } else if (algorithm==ALGORITHM_RKF45) { - ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, true); + ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, true); } else if (algorithm==ALGORITHM_RKDP54) { // only compute k1 on the first step ns_step_rkdp54(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, time==starting_time); } else if (algorithm==ALGORITHM_RKBS32) { // only compute k1 on the first step - ns_step_rkbs32(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, tmp2, tmp3, &tmp4, tmp5, irreversible, time==starting_time); + ns_step_rkbs32(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, &step, &next_step, L, g, fft1, fft2, ifft, &tmp1, tmp2, tmp3, &tmp4, tmp5, irreversible, time==starting_time); } else { fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm); } @@ -649,6 +649,8 @@ int ns_step_rkf45( _Complex double* u, double tolerance, double factor, + double max_delta, + unsigned int adaptive_norm, int K1, int K2, int N1, @@ -722,17 +724,59 @@ int ns_step_rkf45( // compute error err=0; - relative=0; - 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)*(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)])/pow(kx*kx+ky*ky,1.5); + if(adaptive_norm==NORM_L1){ + relative=0; + for(kx=0;kx<=K1;kx++){ + for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){ + 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)*(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)]); + } } } + else if(adaptive_norm==NORM_k3){ + relative=0; + for(kx=0;kx<=K1;kx++){ + for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){ + 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)])/pow(kx*kx+ky*ky,1.5); + } + } + } + else if(adaptive_norm==NORM_k32){ + relative=0; + for(kx=0;kx<=K1;kx++){ + for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){ + 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,0.75); + // 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)])/pow(kx*kx+ky*ky,0.75); + } + } + } + else if(adaptive_norm==NORM_ENSTROPHY){ + double sumu, sumU; + sumu=0; + sumU=0; + for(kx=0;kx<=K1;kx++){ + for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){ + 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)]+tmp[klookup_sym(kx,ky,K2)]); + } + } + err=sqrt(err); + relative=pow((sqrt(sumu)+sqrt(sumU))/sumu, 1./3); + } + else{ + fprintf(stderr,"bug: unknown norm: %u, contact ian.jauslin@rutgers.edu\n", adaptive_norm); + exit(-1); + } // compare relative error with tolerance if(err0 ? -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)]))/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); + if(adaptive_norm==NORM_L1){ + relative=0; + for(kx=0;kx<=K1;kx++){ + for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){ + 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)]); + } } } + else if(adaptive_norm==NORM_k3){ + relative=0; + for(kx=0;kx<=K1;kx++){ + for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){ + 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,1.5); + relative+=cabs(tmp[klookup_sym(kx,ky,K2)]-u[klookup_sym(kx,ky,K2)])/pow(kx*kx+ky*ky,1.5); + } + } + } + else if(adaptive_norm==NORM_k32){ + relative=0; + for(kx=0;kx<=K1;kx++){ + for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){ + 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); + } + } + } + else if(adaptive_norm==NORM_ENSTROPHY){ + double sumu, sumU; + sumu=0; + sumU=0; + 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(tmp[klookup_sym(kx,ky,K2)]); + } + } + err=sqrt(err); + relative=pow((sqrt(sumu)+sqrt(sumU))/sumu, 1./3); + } + else{ + fprintf(stderr,"bug: unknown norm: %u, contact ian.jauslin@rutgers,edu\n", adaptive_norm); + exit(-1); + } // compare relative error with tolerance if(err0 ? -K2 : 1);ky<=K2;ky++){ - // difference between 5th order and 4th order - // 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)])); sumU+=(kx*kx+ky*ky)*cabs2(u[klookup_sym(kx,ky,K2)]+(*delta)*(5179./57600*(*k1)[klookup_sym(kx,ky,K2)]+7571./16695*k3[klookup_sym(kx,ky,K2)]+393./640*k4[klookup_sym(kx,ky,K2)]-92097./339200*k5[klookup_sym(kx,ky,K2)]+187./2100*k6[klookup_sym(kx,ky,K2)]+1./40*(*k2)[klookup_sym(kx,ky,K2)])); sumu+=(kx*kx+ky*ky)*cabs2(tmp[klookup_sym(kx,ky,K2)]); diff --git a/src/navier-stokes.h b/src/navier-stokes.h index 18b6cab..a2b0469 100644 --- a/src/navier-stokes.h +++ b/src/navier-stokes.h @@ -56,11 +56,11 @@ int ns_step_rk4( _Complex double* u, int K1, int K2, int N1, int N2, double nu, // RK2 int ns_step_rk2( _Complex double* u, int K1, int K2, int N1, int N2, double nu, double delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2,fft_vect ifft, _Complex double* tmp1, _Complex double *tmp2, bool irreversible); // Runge-Kutta-Fehlberg -int ns_step_rkf45( _Complex double* u, double tolerance, double factor, int K1, int K2, int N1, int N2, double nu, double* delta, double* next_delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double* k1, _Complex double* k2, _Complex double* k3, _Complex double* k4, _Complex double* k5, _Complex double* k6, _Complex double* tmp, bool irreversible, bool compute_k1); +int ns_step_rkf45( _Complex double* u, double tolerance, double factor, double max_delta, unsigned int adaptive_norm, int K1, int K2, int N1, int N2, double nu, double* delta, double* next_delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double* k1, _Complex double* k2, _Complex double* k3, _Complex double* k4, _Complex double* k5, _Complex double* k6, _Complex double* tmp, bool irreversible, bool compute_k1); // Runge-Kutta-Dromand-Prince int ns_step_rkdp54( _Complex double* u, double tolerance, double factor, double max_delta, unsigned int adaptive_norm, int K1, int K2, int N1, int N2, double nu, double* delta, double* next_delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double** k1, _Complex double** k2, _Complex double* k3, _Complex double* k4, _Complex double* k5, _Complex double* k6, _Complex double* tmp, bool irreversible, bool compute_k1); // Runge-Kutta-Bogacki-Shampine -int ns_step_rkbs32( _Complex double* u, double tolerance, double factor, int K1, int K2, int N1, int N2, double nu, double* delta, double* next_delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double** k1, _Complex double* k2, _Complex double* k3, _Complex double** k4, _Complex double* tmp, bool irreversible, bool compute_k1); +int ns_step_rkbs32( _Complex double* u, double tolerance, double factor, double max_delta, unsigned int adaptive_norm, int K1, int K2, int N1, int N2, double nu, double* delta, double* next_delta, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double** k1, _Complex double* k2, _Complex double* k3, _Complex double** k4, _Complex double* tmp, bool irreversible, bool compute_k1); // right side of Irreversible/reversible Navier-Stokes equation int ns_rhs( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2, double nu, double L, _Complex double* g, fft_vect fft1, fft_vect fft2, fft_vect ifft, bool irreversible);