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

Adaptive timestep: do nto change delta before all the printing

Browse Source
This commit is contained in:
Ian Jauslin 2023-05-16 17:03:58 -04:00
parent cf8cdfe12c
commit dacd5c3651
2 changed files with 54 additions and 33 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

85
src/navier-stokes.c
View File

@ -77,20 +77,25 @@ int uk(
// add 0.1 to ensure proper rounding // add 0.1 to ensure proper rounding
uint64_t n=(uint64_t)((starting_time-fmod(starting_time, print_freq))/print_freq+0.1); uint64_t n=(uint64_t)((starting_time-fmod(starting_time, print_freq))/print_freq+0.1);
// store step (useful for adaptive step methods
double step=delta;
double next_step=step;
// iterate // iterate
time=starting_time; time=starting_time;
while(final_time<0 || time<final_time){ while(final_time<0 || time<final_time){
if(algorithm==ALGORITHM_RK2){ if(algorithm==ALGORITHM_RK2){
ns_step_rk2(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible); ns_step_rk2(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible);
} else if (algorithm==ALGORITHM_RK4) { } else if (algorithm==ALGORITHM_RK4) {
ns_step_rk4(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible); ns_step_rk4(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
} else if (algorithm==ALGORITHM_RKF45) { } else if (algorithm==ALGORITHM_RKF45) {
delta=ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible); 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);
} else { } else {
fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm); fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
} }
time+=delta; time+=step;
step=next_step;
if(time>(n+1)*print_freq){ if(time>(n+1)*print_freq){
n++; n++;
@ -165,6 +170,9 @@ int enstrophy(
// copy initial condition // copy initial condition
copy_u(u, u0, K1, K2); copy_u(u, u0, K1, K2);
// store step (useful for adaptive step methods
double step=delta;
double next_step=step;
// init running average // init running average
avg_a=0; avg_a=0;
@ -180,24 +188,24 @@ int enstrophy(
time=starting_time; time=starting_time;
while(final_time<0 || time<final_time){ while(final_time<0 || time<final_time){
if(algorithm==ALGORITHM_RK2){ if(algorithm==ALGORITHM_RK2){
ns_step_rk2(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible); ns_step_rk2(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible);
} else if (algorithm==ALGORITHM_RK4) { } else if (algorithm==ALGORITHM_RK4) {
ns_step_rk4(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible); ns_step_rk4(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
} else if (algorithm==ALGORITHM_RKF45) { } else if (algorithm==ALGORITHM_RKF45) {
delta=ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible); 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);
} else { } else {
fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm); fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
} }
time+=delta; time+=step;
alpha=compute_alpha(u, K1, K2, g, L); alpha=compute_alpha(u, K1, K2, g, L);
enstrophy=compute_enstrophy(u, K1, K2, L); enstrophy=compute_enstrophy(u, K1, K2, L);
// add to running averages (estimate the total duration of interval as print_freq, will be adjusted later) // add to running averages (estimate the total duration of interval as print_freq, will be adjusted later)
avg_a+=alpha*(delta/print_freq); avg_a+=alpha*(step/print_freq);
avg_en+=enstrophy*(delta/print_freq); avg_en+=enstrophy*(step/print_freq);
avg_en_x_a+=enstrophy*alpha*(delta/print_freq); avg_en_x_a+=enstrophy*alpha*(step/print_freq);
if(time>(n+1)*print_freq){ if(time>(n+1)*print_freq){
n++; n++;
@ -209,8 +217,8 @@ int enstrophy(
// print to stderr so user can follow along // print to stderr so user can follow along
if(algorithm==ALGORITHM_RKF45){ if(algorithm==ALGORITHM_RKF45){
fprintf(stderr,"% .8e % .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",time, avg_a, avg_en, avg_en_x_a, alpha, enstrophy, alpha*enstrophy, delta); fprintf(stderr,"% .8e % .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",time, avg_a, avg_en, avg_en_x_a, alpha, enstrophy, alpha*enstrophy, step);
printf("% .15e % .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",time, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy, delta); printf("% .15e % .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",time, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy, step);
} else { } else {
fprintf(stderr,"% .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",time, avg_a, avg_en, avg_en_x_a, alpha, enstrophy, alpha*enstrophy); fprintf(stderr,"% .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",time, avg_a, avg_en, avg_en_x_a, alpha, enstrophy, alpha*enstrophy);
printf("% .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",time, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy); printf("% .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",time, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy);
@ -224,6 +232,10 @@ int enstrophy(
avg_en_x_a=0; avg_en_x_a=0;
prevtime=time; prevtime=time;
// get ready for next step
step=next_step;
// catch abort signal // catch abort signal
if (g_abort){ if (g_abort){
// print u to stderr if no savefile // print u to stderr if no savefile
@ -254,7 +266,7 @@ int enstrophy(
free(params); free(params);
// write delta if adaptive, and not writing binary // write delta if adaptive, and not writing binary
if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD && (savefile==stderr || savefile==stdout)){ if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD && (savefile==stderr || savefile==stdout)){
fprintf(savefile,";delta=%.15e", delta); fprintf(savefile,";delta=%.15e", step);
} }
// write starting_time if not writing binary // write starting_time if not writing binary
if(savefile==stderr || savefile==stdout){ if(savefile==stderr || savefile==stdout){
@ -273,7 +285,7 @@ int enstrophy(
fwrite(&time, sizeof(double), 1, savefile); fwrite(&time, sizeof(double), 1, savefile);
// extra binary data for adaptive algorithm // extra binary data for adaptive algorithm
if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){ if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){
fwrite(&delta, sizeof(double), 1, savefile); fwrite(&step, sizeof(double), 1, savefile);
} }
} }
} }
@ -319,20 +331,26 @@ int quiet(
// copy initial condition // copy initial condition
copy_u(u, u0, K1, K2); copy_u(u, u0, K1, K2);
// store delta (useful for adaptive step algorithms)
double step=delta;
double next_step=step;
// iterate // iterate
time=starting_time; time=starting_time;
while(final_time<0 || time<final_time){ while(final_time<0 || time<final_time){
if(algorithm==ALGORITHM_RK2){ if(algorithm==ALGORITHM_RK2){
ns_step_rk2(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible); ns_step_rk2(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible);
} else if (algorithm==ALGORITHM_RK4) { } else if (algorithm==ALGORITHM_RK4) {
ns_step_rk4(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible); ns_step_rk4(u, K1, K2, N1, N2, nu, step, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
} else if (algorithm==ALGORITHM_RKF45) { } else if (algorithm==ALGORITHM_RKF45) {
delta=ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible); 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);
} else { } else {
fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm); fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
} }
time+=delta; time+=step;
step=next_step;
} }
// save final entry to savefile // save final entry to savefile
@ -590,7 +608,7 @@ int ns_step_rk2(
// next time step // next time step
// adaptive RK algorithm (Runge-Kutta-Fehlberg) // adaptive RK algorithm (Runge-Kutta-Fehlberg)
double ns_step_rkf45( int ns_step_rkf45(
_Complex double* u, _Complex double* u,
double tolerance, double tolerance,
double factor, double factor,
@ -599,7 +617,8 @@ double ns_step_rkf45(
int N1, int N1,
int N2, int N2,
double nu, double nu,
double delta, double* delta,
double* next_delta,
double L, double L,
_Complex double* g, _Complex double* g,
fft_vect fft1, fft_vect fft1,
@ -623,7 +642,7 @@ double ns_step_rkf45(
// k2 : u(t+1/4*delta) // k2 : u(t+1/4*delta)
for(kx=0;kx<=K1;kx++){ for(kx=0;kx<=K1;kx++){
for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){ for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta/4*k1[klookup_sym(kx,ky,K2)]; tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)/4*k1[klookup_sym(kx,ky,K2)];
} }
} }
ns_rhs(k2, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible); ns_rhs(k2, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
@ -631,7 +650,7 @@ double ns_step_rkf45(
// k3 : u(t+3/8*delta) // k3 : u(t+3/8*delta)
for(kx=0;kx<=K1;kx++){ for(kx=0;kx<=K1;kx++){
for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){ for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta*(3./32*k1[klookup_sym(kx,ky,K2)]+9./32*k2[klookup_sym(kx,ky,K2)]); tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(3./32*k1[klookup_sym(kx,ky,K2)]+9./32*k2[klookup_sym(kx,ky,K2)]);
} }
} }
ns_rhs(k3, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible); ns_rhs(k3, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
@ -639,7 +658,7 @@ double ns_step_rkf45(
// k4 : u(t+12./13*delta) // k4 : u(t+12./13*delta)
for(kx=0;kx<=K1;kx++){ for(kx=0;kx<=K1;kx++){
for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){ for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta*(1932./2197*k1[klookup_sym(kx,ky,K2)]-7200./2197*k2[klookup_sym(kx,ky,K2)]+7296./2197*k3[klookup_sym(kx,ky,K2)]); tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(1932./2197*k1[klookup_sym(kx,ky,K2)]-7200./2197*k2[klookup_sym(kx,ky,K2)]+7296./2197*k3[klookup_sym(kx,ky,K2)]);
} }
} }
ns_rhs(k4, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible); ns_rhs(k4, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
@ -647,7 +666,7 @@ double ns_step_rkf45(
// k5 : u(t+1*delta) // k5 : u(t+1*delta)
for(kx=0;kx<=K1;kx++){ for(kx=0;kx<=K1;kx++){
for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){ for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta*(439./216*k1[klookup_sym(kx,ky,K2)]-8*k2[klookup_sym(kx,ky,K2)]+3680./513*k3[klookup_sym(kx,ky,K2)]-845./4104*k4[klookup_sym(kx,ky,K2)]); tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(439./216*k1[klookup_sym(kx,ky,K2)]-8*k2[klookup_sym(kx,ky,K2)]+3680./513*k3[klookup_sym(kx,ky,K2)]-845./4104*k4[klookup_sym(kx,ky,K2)]);
} }
} }
ns_rhs(k5, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible); ns_rhs(k5, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
@ -655,7 +674,7 @@ double ns_step_rkf45(
// k6 : u(t+1./2*delta) // k6 : u(t+1./2*delta)
for(kx=0;kx<=K1;kx++){ for(kx=0;kx<=K1;kx++){
for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){ for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+delta*(-8./27*k1[klookup_sym(kx,ky,K2)]+2*k2[klookup_sym(kx,ky,K2)]-3544./2565*k3[klookup_sym(kx,ky,K2)]+1859./4104*k4[klookup_sym(kx,ky,K2)]-11./40*k5[klookup_sym(kx,ky,K2)]); tmp[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*delta)*(-8./27*k1[klookup_sym(kx,ky,K2)]+2*k2[klookup_sym(kx,ky,K2)]-3544./2565*k3[klookup_sym(kx,ky,K2)]+1859./4104*k4[klookup_sym(kx,ky,K2)]-11./40*k5[klookup_sym(kx,ky,K2)]);
} }
} }
ns_rhs(k6, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible); ns_rhs(k6, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
@ -666,15 +685,13 @@ double ns_step_rkf45(
for(kx=0;kx<=K1;kx++){ for(kx=0;kx<=K1;kx++){
for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){ for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
// difference between 5th order and 4th order // 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)])); 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 // next step
tmp[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)]); tmp[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)]);
relative+=cabs(tmp[klookup_sym(kx,ky,K2)]); relative+=cabs(tmp[klookup_sym(kx,ky,K2)]);
} }
} }
// new delta
delta=factor*delta*pow(relative*tolerance/err,0.2);
// compare relative error with tolerance // compare relative error with tolerance
if(err<relative*tolerance){ if(err<relative*tolerance){
// add to output // add to output
@ -683,12 +700,16 @@ double ns_step_rkf45(
u[klookup_sym(kx,ky,K2)]+=tmp[klookup_sym(kx,ky,K2)]; u[klookup_sym(kx,ky,K2)]+=tmp[klookup_sym(kx,ky,K2)];
} }
} }
return delta; // next delta to use in future steps
*next_delta=factor*(*delta)*pow(relative*tolerance/err,0.2);
} }
// error too big: repeat with smaller step // error too big: repeat with smaller step
else{ else{
return(ns_step_rkf45(u,tolerance,factor,K1,K2,N1,N2,nu,delta,L,g,fft1,fft2,ifft,k1,k2,k3,k4,k5,k6,tmp,irreversible)); *delta=factor*(*delta)*pow(relative*tolerance/err,0.2);
ns_step_rkf45(u,tolerance,factor,K1,K2,N1,N2,nu,delta,next_delta,L,g,fft1,fft2,ifft,k1,k2,k3,k4,k5,k6,tmp,irreversible);
} }
return 0;
} }
// right side of Irreversible/Reversible Navier-Stokes equation // right side of Irreversible/Reversible Navier-Stokes equation

2
src/navier-stokes.h
View File

@ -56,7 +56,7 @@ int ns_step_rk4( _Complex double* u, int K1, int K2, int N1, int N2, double nu,
// RK2 // 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); 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);
// adaptive RK algorithm (Runge-Kutta-Fehlberg) // adaptive RK algorithm (Runge-Kutta-Fehlberg)
double ns_step_rkf45( _Complex double* u, double tolerance, double factor, 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* k1, _Complex double* k2, _Complex double* k3, _Complex double* k4, _Complex double* k5, _Complex double* k6, _Complex double* tmp, bool irreversible); 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);
// right side of Irreversible/reversible Navier-Stokes equation // 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); 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);