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

keep_en_cst

Browse Source
This commit is contained in:
Ian Jauslin
2023-06-29 16:58:33 -04:00
parent 9c3ad2c400
commit 60e04ac681
4 changed files with 111 additions and 31 deletions
Download Patch File Download Diff File Expand all files Collapse all files

106
src/navier-stokes.c
View File

@@ -40,6 +40,8 @@ int uk(
_Complex double* u0,
_Complex double* g,
bool irreversible,
bool keep_en_cst,
double target_en,
unsigned int algorithm,
double print_freq,
double starting_time,
@@ -88,17 +90,17 @@ int uk(
time=starting_time;
while(final_time<0 || time<final_time){
if(algorithm==ALGORITHM_RK2){
ns_step_rk2(u, K1, K2, N1, N2, nu, step, 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, keep_en_cst, target_en);
} else if (algorithm==ALGORITHM_RK4) {
ns_step_rk4(u, K1, K2, N1, N2, nu, step, 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, keep_en_cst, target_en);
} else if (algorithm==ALGORITHM_RKF45) {
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);
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, keep_en_cst, target_en, 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);
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, keep_en_cst, target_en, time==starting_time);
} else if (algorithm==ALGORITHM_RKBS32) {
// only compute k1 on the first step
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);
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, keep_en_cst, target_en, time==starting_time);
} else {
fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
}
@@ -151,6 +153,8 @@ int enstrophy(
_Complex double* u0,
_Complex double* g,
bool irreversible,
bool keep_en_cst,
double target_en,
unsigned int algorithm,
double print_freq,
double starting_time,
@@ -202,17 +206,17 @@ int enstrophy(
time=starting_time;
while(final_time<0 || time<final_time){
if(algorithm==ALGORITHM_RK2){
ns_step_rk2(u, K1, K2, N1, N2, nu, step, 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, keep_en_cst, target_en);
} else if (algorithm==ALGORITHM_RK4) {
ns_step_rk4(u, K1, K2, N1, N2, nu, step, 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, keep_en_cst, target_en);
} else if (algorithm==ALGORITHM_RKF45) {
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);
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, keep_en_cst, target_en, 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);
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, keep_en_cst, target_en, time==starting_time);
} else if (algorithm==ALGORITHM_RKBS32) {
// only compute k1 on the first step
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);
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, keep_en_cst, target_en, time==starting_time);
} else {
fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
}
@@ -345,6 +349,8 @@ int quiet(
_Complex double* u0,
_Complex double* g,
bool irreversible,
bool keep_en_cst,
double target_en,
unsigned int algorithm,
unsigned int nthreads,
FILE* savefile
@@ -375,17 +381,17 @@ int quiet(
time=starting_time;
while(final_time<0 || time<final_time){
if(algorithm==ALGORITHM_RK2){
ns_step_rk2(u, K1, K2, N1, N2, nu, step, 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, keep_en_cst, target_en);
} else if (algorithm==ALGORITHM_RK4) {
ns_step_rk4(u, K1, K2, N1, N2, nu, step, 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, keep_en_cst, target_en);
} else if (algorithm==ALGORITHM_RKF45) {
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);
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, keep_en_cst, target_en, 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);
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, keep_en_cst, target_en, time==starting_time);
} else if (algorithm==ALGORITHM_RKBS32) {
// only compute k1 on the first step
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);
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, keep_en_cst, target_en, time==starting_time);
} else {
fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers,edu\n",algorithm);
}
@@ -558,7 +564,9 @@ int ns_step_rk4(
_Complex double* tmp1,
_Complex double* tmp2,
_Complex double* tmp3,
bool irreversible
bool irreversible,
bool keep_en_cst,
double target_en
){
int kx,ky;
@@ -616,6 +624,16 @@ int ns_step_rk4(
}
}
// keep enstrophy constant
if(keep_en_cst){
double en=compute_enstrophy(u, K1, K2, L);
for(kx=0;kx<=K1;kx++){
for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
u[klookup_sym(kx,ky,K2)]*=sqrt(target_en/en);
}
}
}
return(0);
}
@@ -635,7 +653,9 @@ int ns_step_rk2(
fft_vect ifft,
_Complex double* tmp1,
_Complex double* tmp2,
bool irreversible
bool irreversible,
bool keep_en_cst,
double target_en
){
int kx,ky;
@@ -657,6 +677,16 @@ int ns_step_rk2(
}
}
// keep enstrophy constant
if(keep_en_cst){
double en=compute_enstrophy(u, K1, K2, L);
for(kx=0;kx<=K1;kx++){
for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
u[klookup_sym(kx,ky,K2)]*=sqrt(target_en/en);
}
}
}
return(0);
}
@@ -688,6 +718,8 @@ int ns_step_rkf45(
_Complex double* k6,
_Complex double* tmp,
bool irreversible,
bool keep_en_cst,
double target_en,
// whether to compute k1 or leave it as is
bool compute_k1
){
@@ -805,12 +837,22 @@ int ns_step_rkf45(
}
// next delta to use in future steps (do not exceed max_delta)
*next_delta=fmin(max_delta, (*delta)*pow(relative*tolerance/err,0.2));
// keep enstrophy constant
if(keep_en_cst){
double en=compute_enstrophy(u, K1, K2, L);
for(kx=0;kx<=K1;kx++){
for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
u[klookup_sym(kx,ky,K2)]*=sqrt(target_en/en);
}
}
}
}
// error too big: repeat with smaller step
else{
*delta=factor*(*delta)*pow(relative*tolerance/err,0.2);
// do not recompute k1
ns_step_rkf45(u,tolerance,factor,max_delta,adaptive_norm,K1,K2,N1,N2,nu,delta,next_delta,L,g,fft1,fft2,ifft,k1,k2,k3,k4,k5,k6,tmp,irreversible,false);
ns_step_rkf45(u,tolerance,factor,max_delta,adaptive_norm,K1,K2,N1,N2,nu,delta,next_delta,L,g,fft1,fft2,ifft,k1,k2,k3,k4,k5,k6,tmp,irreversible,keep_en_cst,target_en,false);
}
return 0;
@@ -844,6 +886,8 @@ int ns_step_rkbs32(
_Complex double** k4,
_Complex double* tmp,
bool irreversible,
bool keep_en_cst,
double target_en,
// whether to compute k1
bool compute_k1
){
@@ -944,12 +988,22 @@ int ns_step_rkbs32(
tmp=*k1;
*k1=*k4;
*k4=tmp;
// keep enstrophy constant
if(keep_en_cst){
double en=compute_enstrophy(u, K1, K2, L);
for(kx=0;kx<=K1;kx++){
for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
u[klookup_sym(kx,ky,K2)]*=sqrt(target_en/en);
}
}
}
}
// error too big: repeat with smaller step
else{
*delta=factor*(*delta)*pow(relative*tolerance/err,1./3);
// this will reuse the same k1 without re-computing it
ns_step_rkbs32(u,tolerance,factor,max_delta,adaptive_norm,K1,K2,N1,N2,nu,delta,next_delta,L,g,fft1,fft2,ifft,k1,k2,k3,k4,tmp,irreversible,false);
ns_step_rkbs32(u,tolerance,factor,max_delta,adaptive_norm,K1,K2,N1,N2,nu,delta,next_delta,L,g,fft1,fft2,ifft,k1,k2,k3,k4,tmp,irreversible,keep_en_cst,target_en,false);
}
return 0;
@@ -984,6 +1038,8 @@ int ns_step_rkdp54(
_Complex double* k6,
_Complex double* tmp,
bool irreversible,
bool keep_en_cst,
double target_en,
// whether to compute k1
bool compute_k1
){
@@ -1109,12 +1165,22 @@ int ns_step_rkdp54(
tmp=*k1;
*k1=*k2;
*k2=tmp;
// keep enstrophy constant
if(keep_en_cst){
double en=compute_enstrophy(u, K1, K2, L);
for(kx=0;kx<=K1;kx++){
for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
u[klookup_sym(kx,ky,K2)]*=sqrt(target_en/en);
}
}
}
}
// error too big: repeat with smaller step
else{
*delta=factor*(*delta)*pow(relative*tolerance/err,1./5);
// this will reuse the same k1 without re-computing it
ns_step_rkdp54(u,tolerance,factor,max_delta,adaptive_norm,K1,K2,N1,N2,nu,delta,next_delta,L,g,fft1,fft2,ifft,k1,k2,k3,k4,k5,k6,tmp,irreversible,false);
ns_step_rkdp54(u,tolerance,factor,max_delta,adaptive_norm,K1,K2,N1,N2,nu,delta,next_delta,L,g,fft1,fft2,ifft,k1,k2,k3,k4,k5,k6,tmp,irreversible,keep_en_cst,target_en,false);
}
return 0;