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

RKF45 algorithm

Browse Source
This commit is contained in:
Ian Jauslin
2023-05-15 20:29:06 -04:00
parent 315e9a98de
commit 5db3680b31
5 changed files with 325 additions and 62 deletions
Download Patch File Download Diff File Expand all files Collapse all files

271
src/navier-stokes.c
View File

@@ -32,26 +32,34 @@ int uk(
double nu,
double delta,
double L,
double adaptive_tolerance,
double adaptive_factor,
_Complex double* u0,
_Complex double* g,
bool irreversible,
unsigned int algorithm,
uint64_t print_freq,
uint64_t starting_time,
uint64_t starting_step,
double starting_time,
unsigned int nthreads,
FILE* savefile
){
double time=starting_time;
_Complex double* u;
_Complex double* tmp1;
_Complex double* tmp2;
_Complex double* tmp3;
_Complex double* tmp4;
_Complex double* tmp5;
_Complex double* tmp6;
_Complex double* tmp7;
uint64_t t;
fft_vect fft1;
fft_vect fft2;
fft_vect ifft;
int kx,ky;
ns_init_tmps(&u, &tmp1, &tmp2, &tmp3, &fft1, &fft2, &ifft, K1, K2, N1, N2, nthreads);
ns_init_tmps(&u, &tmp1, &tmp2, &tmp3, &tmp4, &tmp5, &tmp6, &tmp7, &fft1, &fft2, &ifft, K1, K2, N1, N2, nthreads, algorithm);
// copy initial condition
copy_u(u, u0, K1, K2);
@@ -68,16 +76,22 @@ int uk(
}
// iterate
for(t=starting_time;nsteps==0 || t<starting_time+nsteps;t++){
for(t=starting_step;nsteps==0 || t<starting_step+nsteps;t++){
if(algorithm==ALGORITHM_RK2){
ns_step_rk2(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible);
} else {
} else if (algorithm==ALGORITHM_RK4) {
ns_step_rk4(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
} 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);
} else {
fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
}
time+=delta;
if(t%print_freq==0){
fprintf(stderr,"%lu % .8e ",t,t*delta);
printf("%8lu % .15e ",t,t*delta);
fprintf(stderr,"%lu % .8e ",t,time);
printf("%8lu % .15e ",t,time);
for(kx=-K1;kx<=K1;kx++){
for (ky=-K2;ky<=K2;ky++){
@@ -96,7 +110,7 @@ int uk(
// save final entry to savefile
write_vec_bin(u, K1, K2, savefile);
ns_free_tmps(u, tmp1, tmp2, tmp3, fft1, fft2, ifft);
ns_free_tmps(u, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, fft1, fft2, ifft, algorithm);
return(0);
}
@@ -110,12 +124,15 @@ int enstrophy(
double nu,
double delta,
double L,
double adaptive_tolerance,
double adaptive_factor,
_Complex double* u0,
_Complex double* g,
bool irreversible,
unsigned int algorithm,
uint64_t print_freq,
uint64_t starting_time,
uint64_t starting_step,
double starting_time,
unsigned int nthreads,
FILE* savefile,
// for interrupt recovery
@@ -127,6 +144,11 @@ int enstrophy(
_Complex double* tmp1;
_Complex double* tmp2;
_Complex double* tmp3;
_Complex double* tmp4;
_Complex double* tmp5;
_Complex double* tmp6;
_Complex double* tmp7;
double time=starting_time;
double alpha, enstrophy;
double avg_a,avg_en,avg_en_x_a;
// index
@@ -135,7 +157,7 @@ int enstrophy(
fft_vect fft2;
fft_vect ifft;
ns_init_tmps(&u, &tmp1, &tmp2, &tmp3, &fft1, &fft2, &ifft, K1, K2, N1, N2, nthreads);
ns_init_tmps(&u, &tmp1, &tmp2, &tmp3, &tmp4, &tmp5, &tmp6, &tmp7, &fft1, &fft2, &ifft, K1, K2, N1, N2, nthreads, algorithm);
// copy initial condition
copy_u(u, u0, K1, K2);
@@ -146,26 +168,38 @@ int enstrophy(
avg_en_x_a=0;
// special first case when starting_time is not a multiple of print_freq
uint64_t first_box = print_freq - (starting_time % print_freq);
uint64_t first_box = print_freq - (starting_step % print_freq);
// iterate
for(t=starting_time;nsteps==0 || t<starting_time+nsteps;t++){
for(t=starting_step;nsteps==0 || t<starting_step+nsteps;t++){
if(algorithm==ALGORITHM_RK2){
ns_step_rk2(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible);
} else {
} else if (algorithm==ALGORITHM_RK4) {
ns_step_rk4(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
} 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);
} else {
fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
}
time+=delta;
alpha=compute_alpha(u, K1, K2, g, L);
enstrophy=compute_enstrophy(u, K1, K2, L);
avg_a=average_step(alpha, avg_a, t, starting_time, print_freq, first_box);
avg_en=average_step(enstrophy, avg_en, t, starting_time, print_freq, first_box);
avg_en_x_a=average_step(enstrophy*alpha, avg_en_x_a, t, starting_time, print_freq, first_box);
avg_a=average_step(alpha, avg_a, t, starting_step, print_freq, first_box);
avg_en=average_step(enstrophy, avg_en, t, starting_step, print_freq, first_box);
avg_en_x_a=average_step(enstrophy*alpha, avg_en_x_a, t, starting_step, print_freq, first_box);
if(t>starting_time && t%print_freq==0){
fprintf(stderr,"%lu % .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",t,t*delta, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy);
printf("%8lu % .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",t,t*delta, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy);
if(t>starting_step && t%print_freq==0){
// print to stderr so user can follow along
if(algorithm==ALGORITHM_RKF45){
fprintf(stderr,"%lu % .8e % .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",t,time, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy, delta);
} else {
fprintf(stderr,"%lu % .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",t,time, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy);
}
// print to stdout
printf("%8lu % .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",t,time, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy);
}
// catch abort signal
@@ -189,11 +223,24 @@ int enstrophy(
if(params_string!=NULL) {
char* params=calloc(sizeof(char), strlen(params_string)+1);
strcpy(params, params_string);
remove_entry(params, "starting_step");
remove_entry(params, "starting_time");
remove_entry(params, "init");
remove_entry(params, "nsteps");
fprintf(savefile," -p \"%s;starting_time=%lu;nsteps=%lu;init=file:%s\"", params, t+1, (nsteps+starting_time < t+1 ? 0 : nsteps+starting_time-t-1), savefile_string);
if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){
remove_entry(params, "delta");
}
fprintf(savefile," -p \"%s;starting_step=%lu;nsteps=%lu;init=file:%s", params, t+1, (nsteps+starting_step < t+1 ? 0 : nsteps+starting_step-t-1), savefile_string);
free(params);
// write delta if adaptive, and not writing binary
if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD && (savefile==stderr || savefile==stdout)){
fprintf(savefile,";delta=%.15e;starting_time=%.15e", delta, starting_time);
}
// write starting_time if not adaptive
if(algorithm<ALGORITHM_ADAPTIVE_THRESHOLD){
fprintf(savefile,";starting_time=%.15e", starting_time);
}
fprintf(savefile,"\"");
}
fprintf(savefile," enstrophy\n");
@@ -202,10 +249,17 @@ int enstrophy(
write_vec(u, K1, K2, savefile);
} else {
write_vec_bin(u, K1, K2, savefile);
// extra binary data for adaptive algorithm
if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){
// first binary entry: delta
fwrite(&delta, sizeof(double), 1, savefile);
// ssecond binary entry: starting time
fwrite(&starting_time, sizeof(double), 1, savefile);
}
}
}
ns_free_tmps(u, tmp1, tmp2, tmp3, fft1, fft2, ifft);
ns_free_tmps(u, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, fft1, fft2, ifft, algorithm);
return(0);
}
@@ -219,7 +273,9 @@ int quiet(
double nu,
double delta,
double L,
uint64_t starting_time,
double adaptive_tolerance,
double adaptive_factor,
uint64_t starting_step,
_Complex double* u0,
_Complex double* g,
bool irreversible,
@@ -231,38 +287,50 @@ int quiet(
_Complex double* tmp1;
_Complex double* tmp2;
_Complex double* tmp3;
_Complex double* tmp4;
_Complex double* tmp5;
_Complex double* tmp6;
_Complex double* tmp7;
uint64_t t;
fft_vect fft1;
fft_vect fft2;
fft_vect ifft;
ns_init_tmps(&u, &tmp1, &tmp2, &tmp3, &fft1, &fft2, &ifft, K1, K2, N1, N2, nthreads);
ns_init_tmps(&u, &tmp1, &tmp2, &tmp3, &tmp4, &tmp5, &tmp6, &tmp7, &fft1, &fft2, &ifft, K1, K2, N1, N2, nthreads, algorithm);
// copy initial condition
copy_u(u, u0, K1, K2);
// iterate
for(t=starting_time;nsteps==0 || t<starting_time+nsteps;t++){
for(t=starting_step;nsteps==0 || t<starting_step+nsteps;t++){
if(algorithm==ALGORITHM_RK2){
ns_step_rk2(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible);
} else {
} else if (algorithm==ALGORITHM_RK4) {
ns_step_rk4(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
} 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);
} else {
fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
}
}
// save final entry to savefile
write_vec(u, K1, K2, savefile);
ns_free_tmps(u, tmp1, tmp2, tmp3, fft1, fft2, ifft);
ns_free_tmps(u, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, fft1, fft2, ifft, algorithm);
return(0);
}
// initialize vectors for computation
int ns_init_tmps(
int ns_init_tmps(
_Complex double ** u,
_Complex double ** tmp1,
_Complex double ** tmp2,
_Complex double ** tmp3,
_Complex double ** tmp4,
_Complex double ** tmp5,
_Complex double ** tmp6,
_Complex double ** tmp7,
fft_vect* fft1,
fft_vect* fft2,
fft_vect* ifft,
@@ -270,15 +338,31 @@ int ns_init_tmps(
int K2,
int N1,
int N2,
unsigned int nthreads
unsigned int nthreads,
unsigned int algorithm
){
// velocity field
*u=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
// allocate tmp vectors for computation
*tmp1=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
*tmp2=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
*tmp3=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
if(algorithm==ALGORITHM_RK2){
*tmp1=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
*tmp2=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
} else if (algorithm==ALGORITHM_RK4){
*tmp1=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
*tmp2=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
*tmp3=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
} else if (algorithm==ALGORITHM_RKF45){
*tmp1=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
*tmp2=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
*tmp3=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
*tmp4=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
*tmp5=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
*tmp6=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
*tmp7=calloc(sizeof(_Complex double),K1*(2*K2+1)+K2);
} else {
fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
};
// init threads
fftw_init_threads();
@@ -301,9 +385,14 @@ int ns_free_tmps(
_Complex double* tmp1,
_Complex double* tmp2,
_Complex double* tmp3,
_Complex double* tmp4,
_Complex double* tmp5,
_Complex double* tmp6,
_Complex double* tmp7,
fft_vect fft1,
fft_vect fft2,
fft_vect ifft
fft_vect ifft,
unsigned int algorithm
){
// free memory
fftw_destroy_plan(fft1.fft_plan);
@@ -315,9 +404,24 @@ int ns_free_tmps(
fftw_cleanup_threads();
free(tmp3);
free(tmp2);
free(tmp1);
if(algorithm==ALGORITHM_RK2){
free(tmp1);
free(tmp2);
} else if (algorithm==ALGORITHM_RK4){
free(tmp1);
free(tmp2);
free(tmp3);
} else if (algorithm==ALGORITHM_RKF45){
free(tmp1);
free(tmp2);
free(tmp3);
free(tmp4);
free(tmp5);
free(tmp6);
free(tmp7);
} else {
fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
};
free(u);
@@ -462,6 +566,105 @@ int ns_step_rk2(
return(0);
}
// next time step
// 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 kx,ky;
double err;
// k1: u(t)
ns_rhs(k1, u, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
// k2 : u(t+1/4*delta)
for(kx=0;kx<=K1;kx++){
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)];
}
}
ns_rhs(k2, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
// k3 : u(t+3/8*delta)
for(kx=0;kx<=K1;kx++){
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)]);
}
}
ns_rhs(k3, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
// k4 : u(t+12./13*delta)
for(kx=0;kx<=K1;kx++){
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)]);
}
}
ns_rhs(k4, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
// k5 : u(t+1*delta)
for(kx=0;kx<=K1;kx++){
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)]);
}
}
ns_rhs(k5, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
// k6 : u(t+1./2*delta)
for(kx=0;kx<=K1;kx++){
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)]);
}
}
ns_rhs(k6, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
// compute error
err=0;
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)]));
}
}
// new delta
delta=factor*delta*pow(tolerance/err,0.2);
// compare relative error with tolerance
if(err<tolerance){
// add to output
for(kx=0;kx<=K1;kx++){
for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
u[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)]);
}
}
return delta;
}
// error too big: repeat with smaller step
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));
}
}
// right side of Irreversible/Reversible Navier-Stokes equation
int ns_rhs(
_Complex double* out,