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

Rewrite cost function for adaptive step

Browse Source
This commit is contained in:
Ian Jauslin
2024-11-18 17:09:17 -05:00
parent 9fa10c8db4
commit d81ad18618
6 changed files with 211 additions and 264 deletions
Download Patch File Download Diff File Expand all files Collapse all files

9
src/constants.cpp
View File

@ -40,7 +40,8 @@ limitations under the License.
#define ALGORITHM_RKDP54 1002
#define ALGORITHM_RKBS32 1003
#define NORM_L1 1
#define NORM_k3 2
#define NORM_k32 3
#define NORM_ENSTROPHY 4
#define COST_L1 1
#define COST_k3 2
#define COST_k32 3
#define COST_ENSTROPHY 4
#define COST_ALPHA 5

48
src/main.c
View File

@ -46,7 +46,7 @@ typedef struct nstrophy_parameters {
double adaptive_tolerance;
double adaptive_factor;
double max_delta;
unsigned int adaptive_norm;
unsigned int adaptive_cost;
double print_freq;
int seed;
double starting_time;
@ -272,19 +272,19 @@ int main (
// run command
if (command==COMMAND_UK){
uk(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_norm, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_en, parameters.algorithm, parameters.print_freq, parameters.starting_time, nthreads, savefile);
uk(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_cost, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_en, parameters.algorithm, parameters.print_freq, parameters.starting_time, nthreads, savefile);
}
else if(command==COMMAND_ENSTROPHY){
// register signal handler to handle aborts
signal(SIGINT, sig_handler);
signal(SIGTERM, sig_handler);
enstrophy(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_norm, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_en, parameters.algorithm, parameters.print_freq, parameters.starting_time, parameters.print_alpha, nthreads, savefile, utfile, (char*)argv[0], dstring_to_str_noinit(&param_str), dstring_to_str_noinit(&savefile_str), dstring_to_str_noinit(&utfile_str));
enstrophy(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_cost, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_en, parameters.algorithm, parameters.print_freq, parameters.starting_time, parameters.print_alpha, nthreads, savefile, utfile, (char*)argv[0], dstring_to_str_noinit(&param_str), dstring_to_str_noinit(&savefile_str), dstring_to_str_noinit(&utfile_str));
}
else if(command==COMMAND_QUIET){
quiet(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_norm, parameters.starting_time, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_en, parameters.algorithm, nthreads, savefile);
quiet(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_cost, parameters.starting_time, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_en, parameters.algorithm, nthreads, savefile);
}
else if(command==COMMAND_LYAPUNOV){
lyapunov(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.lyapunov_reset, parameters.nu, parameters.D_epsilon, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_norm, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_en, parameters.algorithm, parameters.starting_time, nthreads);
lyapunov(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.lyapunov_reset, parameters.nu, parameters.D_epsilon, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_cost, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_en, parameters.algorithm, parameters.starting_time, nthreads);
}
else if(command==0){
fprintf(stderr, "error: no command specified\n");
@ -393,18 +393,21 @@ int print_params(
}
if(parameters.algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){
switch(parameters.adaptive_norm){
case NORM_L1:
fprintf(file,", norm=L1");
switch(parameters.adaptive_cost){
case COST_L1:
fprintf(file,", cost=L1");
break;
case NORM_k3:
fprintf(file,", norm=k3");
case COST_k3:
fprintf(file,", cost=k3");
break;
case NORM_k32:
fprintf(file,", norm=k32");
case COST_k32:
fprintf(file,", cost=k32");
break;
case NORM_ENSTROPHY:
fprintf(file,", norm=enstrophy");
case COST_ENSTROPHY:
fprintf(file,", cost=enstrophy");
break;
case COST_ALPHA:
fprintf(file,", cost=alpha");
break;
}
}
@ -544,7 +547,7 @@ int set_default_params(
parameters->adaptive_tolerance=1e-11;
parameters->adaptive_factor=0.9;
parameters->max_delta=1e-2;
parameters->adaptive_norm=NORM_L1;
parameters->adaptive_cost=COST_L1;
parameters->final_time=100000;
parameters->print_freq=1;
parameters->starting_time=0;
@ -759,21 +762,24 @@ int set_parameter(
return(-1);
}
}
else if (strcmp(lhs,"adaptive_norm")==0){
else if (strcmp(lhs,"adaptive_cost")==0){
if (strcmp(rhs,"L1")==0){
parameters->adaptive_norm=NORM_L1;
parameters->adaptive_cost=COST_L1;
}
else if (strcmp(rhs,"k3")==0){
parameters->adaptive_norm=NORM_k3;
parameters->adaptive_cost=COST_k3;
}
else if (strcmp(rhs,"k32")==0){
parameters->adaptive_norm=NORM_k32;
parameters->adaptive_cost=COST_k32;
}
else if (strcmp(rhs,"enstrophy")==0){
parameters->adaptive_norm=NORM_ENSTROPHY;
parameters->adaptive_cost=COST_ENSTROPHY;
}
else if (strcmp(rhs,"alpha")==0){
parameters->adaptive_cost=COST_ALPHA;
}
else{
fprintf(stderr, "error: unrecognized adaptive_norm '%s'\n",rhs);
fprintf(stderr, "error: unrecognized adaptive_cost '%s'\n",rhs);
return(-1);
}
}

307
src/navier-stokes.c
View File

@ -36,7 +36,7 @@ int uk(
double adaptive_tolerance,
double adaptive_factor,
double max_delta,
unsigned int adaptive_norm,
unsigned int adaptive_cost,
_Complex double* u0,
_Complex double* g,
bool irreversible,
@ -90,7 +90,7 @@ int uk(
time=starting_time;
while(final_time<0 || time<final_time){
// step
ns_step(algorithm, u, K1, K2, N1, N2, nu, &step, &next_step, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, L, g, time, starting_time, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en);
ns_step(algorithm, u, K1, K2, N1, N2, nu, &step, &next_step, adaptive_tolerance, adaptive_factor, max_delta, adaptive_cost, L, g, time, starting_time, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en);
time+=step;
step=next_step;
@ -136,7 +136,7 @@ int enstrophy(
double adaptive_tolerance,
double adaptive_factor,
double max_delta,
unsigned int adaptive_norm,
unsigned int adaptive_cost,
_Complex double* u0,
_Complex double* g,
bool irreversible,
@ -193,7 +193,7 @@ int enstrophy(
// iterate
time=starting_time;
while(final_time<0 || time<final_time){
ns_step(algorithm, u, K1, K2, N1, N2, nu, &step, &next_step, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, L, g, time, starting_time, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en);
ns_step(algorithm, u, K1, K2, N1, N2, nu, &step, &next_step, adaptive_tolerance, adaptive_factor, max_delta, adaptive_cost, L, g, time, starting_time, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en);
time+=step;
@ -281,7 +281,7 @@ int quiet(
double adaptive_tolerance,
double adaptive_factor,
double max_delta,
unsigned int adaptive_norm,
unsigned int adaptive_cost,
double starting_time,
_Complex double* u0,
_Complex double* g,
@ -317,7 +317,7 @@ int quiet(
// iterate
time=starting_time;
while(final_time<0 || time<final_time){
ns_step(algorithm, u, K1, K2, N1, N2, nu, &step, &next_step, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, L, g, time, starting_time, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en);
ns_step(algorithm, u, K1, K2, N1, N2, nu, &step, &next_step, adaptive_tolerance, adaptive_factor, max_delta, adaptive_cost, L, g, time, starting_time, fft1, fft2, ifft, &tmp1, &tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en);
time+=step;
step=next_step;
@ -483,7 +483,7 @@ int ns_step(
double adaptive_tolerance,
double adaptive_factor,
double max_delta,
unsigned int adaptive_norm,
unsigned int adaptive_cost,
double L,
_Complex double* g,
double time,
@ -508,15 +508,15 @@ int ns_step(
} else if (algorithm==ALGORITHM_RK4) {
ns_step_rk4(u, K1, K2, N1, N2, nu, *delta, 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, delta, next_delta, L, g, fft1, fft2, ifft, *tmp1, *tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en, true);
ns_step_rkf45(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_cost, K1, K2, N1, N2, nu, delta, next_delta, 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
// first-same-as-last with 2-nd argument
ns_step_rkdp54(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, delta, next_delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, irreversible, keep_en_cst, target_en, time==starting_time);
ns_step_rkdp54(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_cost, K1, K2, N1, N2, nu, delta, next_delta, 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
// first-same-as-last with 4-th argument
ns_step_rkbs32(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_norm, K1, K2, N1, N2, nu, delta, next_delta, L, g, fft1, fft2, ifft, tmp1, tmp3, tmp4, tmp2, tmp5, irreversible, keep_en_cst, target_en, time==starting_time);
ns_step_rkbs32(u, adaptive_tolerance, adaptive_factor, max_delta, adaptive_cost, K1, K2, N1, N2, nu, delta, next_delta, L, g, fft1, fft2, ifft, tmp1, tmp3, tmp4, tmp2, tmp5, irreversible, keep_en_cst, target_en, time==starting_time);
} else {
fprintf(stderr,"bug: unknown algorithm: %u, contact ian.jauslin@rutgers.edu\n",algorithm);
}
@ -674,7 +674,7 @@ int ns_step_rkf45(
double tolerance,
double factor,
double max_delta,
unsigned int adaptive_norm,
unsigned int adaptive_cost,
int K1,
int K2,
int N1,
@ -701,7 +701,7 @@ int ns_step_rkf45(
bool compute_k1
){
int kx,ky;
double err,relative;
double cost;
// k1: u(t)
if(compute_k1){
@ -748,72 +748,29 @@ int ns_step_rkf45(
}
ns_rhs(k6, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
// compute error
err=0;
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)]);
}
// next step
for(kx=0;kx<=K1;kx++){
for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
// u
tmp[klookup_sym(kx,ky,K2)]=u[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)]);
// U: save to k6, which is no longer needed
k6[klookup_sym(kx,ky,K2)]=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)]);
}
}
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(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)]));
sumu+=(kx*kx+ky*ky)*cabs2(u[klookup_sym(kx,ky,K2)]+tmp[klookup_sym(kx,ky,K2)]);
}
}
err=sqrt(err);
relative=(sqrt(sumu)+sqrt(sumU))/sumu;
}
else{
fprintf(stderr,"bug: unknown norm: %u, contact ian.jauslin@rutgers.edu\n", adaptive_norm);
exit(-1);
}
// cost function
cost=ns_adaptive_cost(tmp,k6,adaptive_cost,K1,K2,g,L);
// compare relative error with tolerance
if(err<relative*tolerance){
// add to output
if(cost<tolerance){
// copy to output
for(kx=0;kx<=K1;kx++){
for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
u[klookup_sym(kx,ky,K2)]+=tmp[klookup_sym(kx,ky,K2)];
u[klookup_sym(kx,ky,K2)]=tmp[klookup_sym(kx,ky,K2)];
}
}
// next delta to use in future steps (do not exceed max_delta)
*next_delta=fmin(max_delta, (*delta)*pow(relative*tolerance/err,0.2));
*next_delta=fmin(max_delta, (*delta)*pow(tolerance/cost,0.2));
// keep enstrophy constant
if(keep_en_cst){
@ -827,9 +784,9 @@ int ns_step_rkf45(
}
// error too big: repeat with smaller step
else{
*delta=factor*(*delta)*pow(relative*tolerance/err,0.2);
*delta=factor*(*delta)*pow(tolerance/cost,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,keep_en_cst,target_en,false);
ns_step_rkf45(u,tolerance,factor,max_delta,adaptive_cost,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;
@ -843,7 +800,7 @@ int ns_step_rkbs32(
double tolerance,
double factor,
double max_delta,
unsigned int adaptive_norm,
unsigned int adaptive_cost,
int K1,
int K2,
int N1,
@ -869,7 +826,7 @@ int ns_step_rkbs32(
bool compute_k1
){
int kx,ky;
double err,relative;
double cost;
// k1: u(t)
// only compute it if it is the first step (otherwise, it has already been computed due to the First Same As Last property)
@ -894,64 +851,27 @@ int ns_step_rkbs32(
ns_rhs(k3, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
// k4 : u(t+delta)
// tmp cpmputed here is the next step
// tmp computed here is the next step
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)*(2./9*(*k1)[klookup_sym(kx,ky,K2)]+1./3*k2[klookup_sym(kx,ky,K2)]+4./9*k3[klookup_sym(kx,ky,K2)]);
}
}
ns_rhs(*k4, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
// compute error
err=0;
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(u[klookup_sym(kx,ky,K2)]+(*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=(sqrt(sumu)+sqrt(sumU))/sumu;
}
else{
fprintf(stderr,"bug: unknown norm: %u, contact ian.jauslin@rutgers,edu\n", adaptive_norm);
exit(-1);
}
// compare relative error with tolerance
if(err<relative*tolerance){
// next step
for(kx=0;kx<=K1;kx++){
for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
// U: store in k3, which is no longer needed
k3[klookup_sym(kx,ky,K2)]=u[klookup_sym(kx,ky,K2)]+(*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)]);
}
}
// compute cost
cost=ns_adaptive_cost(tmp, k3, adaptive_cost, K1, K2, g, L);
// compare cost with tolerance
if(cost<tolerance){
// add to output
for(kx=0;kx<=K1;kx++){
for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
@ -959,7 +879,7 @@ int ns_step_rkbs32(
}
}
// next delta to use in future steps (do not exceed max_delta)
*next_delta=fmin(max_delta, (*delta)*pow(relative*tolerance/err,1./3));
*next_delta=fmin(max_delta, (*delta)*pow(tolerance/cost,1./3));
// k1 in the next step will be this k4 (first same as last)
tmp=*k1;
@ -978,9 +898,9 @@ int ns_step_rkbs32(
}
// error too big: repeat with smaller step
else{
*delta=factor*(*delta)*pow(relative*tolerance/err,1./3);
*delta=factor*(*delta)*pow(tolerance/cost,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,keep_en_cst,target_en,false);
ns_step_rkbs32(u,tolerance,factor,max_delta,adaptive_cost,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;
@ -993,7 +913,7 @@ int ns_step_rkdp54(
double tolerance,
double factor,
double max_delta,
unsigned int adaptive_norm,
unsigned int adaptive_cost,
int K1,
int K2,
int N1,
@ -1021,7 +941,7 @@ int ns_step_rkdp54(
bool compute_k1
){
int kx,ky;
double err,relative;
double cost;
// k1: u(t)
// only compute it if it is the first step (otherwise, it has already been computed due to the First Same As Last property)
@ -1078,57 +998,20 @@ int ns_step_rkdp54(
}
// store in k2, which is not needed anymore
ns_rhs(*k2, tmp, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft, irreversible);
//next step
for(kx=0;kx<=K1;kx++){
for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
// U: store in k6, which is not needed anymore
k6[klookup_sym(kx,ky,K2)]=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)]);
}
}
// compute error
err=0;
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)*(-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)]);
}
}
}
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)*(-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);
}
}
}
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)*(-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,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)*(-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)]);
}
}
err=sqrt(err);
relative=(sqrt(sumu)+sqrt(sumU))/sumu;
}
else{
fprintf(stderr,"bug: unknown norm: %u, contact ian.jauslin@rutgers,edu\n", adaptive_norm);
exit(-1);
}
// compute cost
cost=ns_adaptive_cost(tmp, k6, adaptive_cost, K1, K2, g, L);
// compare relative error with tolerance
if(err<relative*tolerance){
if(cost<tolerance){
// add to output
for(kx=0;kx<=K1;kx++){
for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
@ -1136,9 +1019,9 @@ int ns_step_rkdp54(
}
}
// next delta to use in future steps (do not exceed max_delta)
*next_delta=fmin(max_delta, (*delta)*pow(relative*tolerance/err,1./5));
*next_delta=fmin(max_delta, (*delta)*pow(tolerance/cost,0.2));
// k1 in the next step will be this k4 (first same as last)
// k1 in the next step will be this k7 (first same as last)
tmp=*k1;
*k1=*k2;
*k2=tmp;
@ -1155,9 +1038,71 @@ int ns_step_rkdp54(
}
// error too big: repeat with smaller step
else{
*delta=factor*(*delta)*pow(relative*tolerance/err,1./5);
*delta=factor*(*delta)*pow(tolerance/cost,0.2);
// 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,keep_en_cst,target_en,false);
ns_step_rkdp54(u,tolerance,factor,max_delta,adaptive_cost,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;
}
// compute error for adaptive step methods
double ns_adaptive_cost(
_Complex double* u,
_Complex double* U,
unsigned int adaptive_cost,
int K1,
int K2,
_Complex double* g,
double L
){
int kx,ky;
if(adaptive_cost==COST_L1){
double cost=0;
double relative=0;
for(kx=0;kx<=K1;kx++){
for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
cost+=cabs(u[klookup_sym(kx,ky,K2)]-U[klookup_sym(kx,ky,K2)]);
relative+=cabs(u[klookup_sym(kx,ky,K2)]);
}
}
return cost/relative;
}
else if(adaptive_cost==COST_k3){
double cost=0;
double relative=0;
for(kx=0;kx<=K1;kx++){
for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
cost+=cabs(u[klookup_sym(kx,ky,K2)]-U[klookup_sym(kx,ky,K2)])/pow(kx*kx+ky*ky,3);
relative+=cabs(u[klookup_sym(kx,ky,K2)])/pow(kx*kx+ky*ky,3);
}
}
return cost/relative;
}
else if(adaptive_cost==COST_k32){
double cost=0;
double relative=0;
for(kx=0;kx<=K1;kx++){
for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
cost+=cabs(u[klookup_sym(kx,ky,K2)]-U[klookup_sym(kx,ky,K2)])/pow(kx*kx+ky*ky,1.5);
relative+=cabs(u[klookup_sym(kx,ky,K2)])/pow(kx*kx+ky*ky,1.5);
}
}
return cost/relative;
}
else if(adaptive_cost==COST_ENSTROPHY){
double enu=compute_enstrophy(u,K1,K2,L);
return fabs(enu-compute_enstrophy(U,K1,K2,L))/enu;
}
else if(adaptive_cost==COST_ALPHA){
double alu=compute_alpha(u,K1,K2,g,L);
return fabs((alu-compute_alpha(U,K1,K2,g,L))/alu);
}
else{
fprintf(stderr,"bug: unknown norm: %u, contact ian.jauslin@rutgers.edu\n", adaptive_cost);
exit(-1);
}
return 0;

18
src/navier-stokes.h
View File

@ -31,13 +31,13 @@ typedef struct fft_vects {
} fft_vect;
// compute u_k
int uk( int K1, int K2, int N1, int N2, double final_time, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_norm, _Complex double* u0, _Complex double* g, bool irreversible, bool keep_en_cst, double target_en, unsigned int algorithm, double print_freq, double starting_time, unsigned int nthreadsl, FILE* savefile);
int uk( int K1, int K2, int N1, int N2, double final_time, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_cost, _Complex double* u0, _Complex double* g, bool irreversible, bool keep_en_cst, double target_en, unsigned int algorithm, double print_freq, double starting_time, unsigned int nthreadsl, FILE* savefile);
// compute enstrophy and alpha
int enstrophy( int K1, int K2, int N1, int N2, double final_time, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_norm, _Complex double* u0, _Complex double* g, bool irreversible, bool keep_en_cst, double target_en, unsigned int algorithm, double print_freq, double starting_time, bool print_alpha, unsigned int nthreads, FILE* savefile, FILE* utfile, const char* cmd_string, const char* params_string, const char* savefile_string, const char* utfile_string);
int enstrophy( int K1, int K2, int N1, int N2, double final_time, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_cost, _Complex double* u0, _Complex double* g, bool irreversible, bool keep_en_cst, double target_en, unsigned int algorithm, double print_freq, double starting_time, bool print_alpha, unsigned int nthreads, FILE* savefile, FILE* utfile, const char* cmd_string, const char* params_string, const char* savefile_string, const char* utfile_string);
// compute solution as a function of time, but do not print anything (useful for debugging)
int quiet( int K1, int K2, int N1, int N2, double final_time, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_norm, double starting_time, _Complex double* u0, _Complex double* g, bool irreversible, bool keep_en_cst, double target_en, unsigned int algorithm, unsigned int nthreads, FILE* savefile);
int quiet( int K1, int K2, int N1, int N2, double final_time, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_cost, double starting_time, _Complex double* u0, _Complex double* g, bool irreversible, bool keep_en_cst, double target_en, unsigned int algorithm, unsigned int nthreads, FILE* savefile);
// initialize vectors for computation
@ -49,17 +49,21 @@ int ns_free_tmps( _Complex double* u, _Complex double* tmp1, _Complex double *tm
int copy_u( _Complex double* u, _Complex double* u0, int K1, int K2);
// next time step for Irreversible/reversible Navier-Stokes equation
int ns_step( unsigned int algorithm, _Complex double* u, int K1, int K2, int N1, int N2, double nu, double* delta, double* next_delta, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_norm, double L, _Complex double* g, double time, double starting_time, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double** tmp1, _Complex double** tmp2, _Complex double* tmp3, _Complex double* tmp4, _Complex double* tmp5, _Complex double* tmp6, _Complex double* tmp7, bool irreversible, bool keep_en_cst, double target_en);
int ns_step( unsigned int algorithm, _Complex double* u, int K1, int K2, int N1, int N2, double nu, double* delta, double* next_delta, double adaptive_tolerance, double adaptive_factor, double max_delta, unsigned int adaptive_cost, double L, _Complex double* g, double time, double starting_time, fft_vect fft1, fft_vect fft2, fft_vect ifft, _Complex double** tmp1, _Complex double** tmp2, _Complex double* tmp3, _Complex double* tmp4, _Complex double* tmp5, _Complex double* tmp6, _Complex double* tmp7, bool irreversible, bool keep_en_cst, double target_en);
// RK4
int ns_step_rk4( _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, _Complex double *tmp3, bool irreversible, bool keep_en_cst, double target_en);
// 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, bool keep_en_cst, double target_en);
// Runge-Kutta-Fehlberg
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 keep_en_cst, double target_en, bool compute_k1);
int ns_step_rkf45( _Complex double* u, double tolerance, double factor, double max_delta, unsigned int adaptive_cost, 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 keep_en_cst, double target_en, 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 keep_en_cst, double target_en, bool compute_k1);
int ns_step_rkdp54( _Complex double* u, double tolerance, double factor, double max_delta, unsigned int adaptive_cost, 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 keep_en_cst, double target_en, bool compute_k1);
// Runge-Kutta-Bogacki-Shampine
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 keep_en_cst, double target_en, bool compute_k1);
int ns_step_rkbs32( _Complex double* u, double tolerance, double factor, double max_delta, unsigned int adaptive_cost, 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 keep_en_cst, double target_en, bool compute_k1);
// cost function for the adaptive iterations
double ns_adaptive_cost( _Complex double* u, _Complex double* U, unsigned int adaptive_cost, int K1, int K2, _Complex double* g, double L);
// 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);