Select norm in parameter

This commit is contained in:
Ian Jauslin 2023-06-13 23:56:35 -04:00
parent 3b58896e3b
commit 251426aaaa
6 changed files with 112 additions and 24 deletions

View File

@ -108,6 +108,11 @@ should be a `;` sperated list of `key=value` pairs. The possible keys are
* `max_delta` (double, default 1e-2): when using the adaptive step, do not
exceet `delta_max`.
* `adaptive_norm`: norm to use to estimate the error of the adaptive method:
`L1` (default) for the normalized L1 norm, `k3` for the normalized L1 norm of
f_k/|k|^3, `k32` for the normalized L1 norm, `enstrophy` for the enstrophy
norm.
# Interrupting/resuming the computation

View File

@ -145,6 +145,7 @@ To be safe, we also set a maximal value for $\delta$ via the {\tt max\_delta} pa
\indent
The choice of the norm $\|\cdot\|$ matters.
It can be made by specifying the parameter {\tt adaptive\_norm}.
\begin{itemize}
\item A naive choice is to take $\|\cdot\|$ to be the normalized $L_1$ norm:
\begin{equation}
@ -154,6 +155,7 @@ The choice of the norm $\|\cdot\|$ matters.
\mathcal N:=\sum_k|\hat u_k^{(n)}-\hat u_k^{(n-1)}|
.
\end{equation}
This norm is selected by choosing {\tt adaptive\_norm=L1}.
\item Empirically, we have found that $|\hat u-\hat U|$ behaves like $k^{-3}$ for {\tt RKDP54} and {\tt RKF45}, and like $k^{-\frac32}$ for {\tt RKBS32}, so a norm of the form
\begin{equation}
@ -168,6 +170,7 @@ The choice of the norm $\|\cdot\|$ matters.
\mathcal N:=\sum_k|\hat u_k^{(n)}-\hat u_k^{(n-1)}|k^{-\frac32}
\end{equation}
are sensible choices.
These norms are selected by choosing {\tt adaptive\_norm=k3} and {\tt adaptive\_norm=k32} respectively.
\item
Another option is to define a norm based on the expression of the enstrophy\-~(\ref{enstrophy}):
@ -193,6 +196,7 @@ The choice of the norm $\|\cdot\|$ matters.
\|\hat u-\hat U\|
.
\end{equation}
This norm is selected by choosing {\tt adaptive\_norm=enstrophy}.
\end{itemize}
\bigskip

View File

@ -37,3 +37,7 @@ 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

View File

@ -48,6 +48,7 @@ typedef struct nstrophy_parameters {
double adaptive_tolerance;
double adaptive_factor;
double max_delta;
unsigned int adaptive_norm;
double print_freq;
int seed;
double starting_time;
@ -171,16 +172,16 @@ 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, u0, g, parameters.irreversible, 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_norm, u0, g, parameters.irreversible, 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, u0, g, parameters.irreversible, parameters.algorithm, parameters.print_freq, parameters.starting_time, nthreads, savefile, (char*)argv[0], param_str, savefile_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_norm, u0, g, parameters.irreversible, parameters.algorithm, parameters.print_freq, parameters.starting_time, nthreads, savefile, (char*)argv[0], param_str, savefile_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.starting_time, u0, g, parameters.irreversible, 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_norm, parameters.starting_time, u0, g, parameters.irreversible, parameters.algorithm, nthreads, savefile);
}
else if(command==0){
fprintf(stderr, "error: no command specified\n");
@ -278,6 +279,23 @@ int print_params(
break;
}
if(parameters.algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){
switch(parameters.adaptive_norm){
case NORM_L1:
fprintf(file,", norm=L1");
break;
case NORM_k3:
fprintf(file,", norm=k3");
break;
case NORM_k32:
fprintf(file,", norm=k32");
break;
case NORM_ENSTROPHY:
fprintf(file,", norm=enstrophy");
break;
}
}
fprintf(file,"\n");
return 0;
@ -399,6 +417,7 @@ int read_params(
parameters->adaptive_tolerance=1e-11;
parameters->adaptive_factor=0.9;
parameters->max_delta=1e-2;
parameters->adaptive_norm=NORM_L1;
parameters->final_time=100000;
parameters->print_freq=1;
parameters->starting_time=0;
@ -606,6 +625,24 @@ int set_parameter(
return(-1);
}
}
else if (strcmp(lhs,"adaptive_norm")==0){
if (strcmp(rhs,"L1")==0){
parameters->adaptive_norm=NORM_L1;
}
else if (strcmp(rhs,"k3")==0){
parameters->adaptive_norm=NORM_k3;
}
else if (strcmp(rhs,"k32")==0){
parameters->adaptive_norm=NORM_k32;
}
else if (strcmp(rhs,"enstrophy")==0){
parameters->adaptive_norm=NORM_ENSTROPHY;
}
else{
fprintf(stderr, "error: unrecognized adaptive_norm '%s'\n",rhs);
return(-1);
}
}
else if (strcmp(lhs,"print_freq")==0){
ret=sscanf(rhs,"%lf",&(parameters->print_freq));
if(ret!=1){

View File

@ -36,6 +36,7 @@ int uk(
double adaptive_tolerance,
double adaptive_factor,
double max_delta,
unsigned int adaptive_norm,
_Complex double* u0,
_Complex double* g,
bool irreversible,
@ -94,7 +95,7 @@ int uk(
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);
} else if (algorithm==ALGORITHM_RKDP54) {
// only compute k1 on the first step
ns_step_rkdp54(u, adaptive_tolerance, adaptive_factor, max_delta, 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, 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);
@ -145,6 +146,7 @@ int enstrophy(
double adaptive_tolerance,
double adaptive_factor,
double max_delta,
unsigned int adaptive_norm,
_Complex double* u0,
_Complex double* g,
bool irreversible,
@ -204,7 +206,7 @@ int enstrophy(
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);
} else if (algorithm==ALGORITHM_RKDP54) {
// only compute k1 on the first step
ns_step_rkdp54(u, adaptive_tolerance, adaptive_factor, max_delta, 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, 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);
@ -318,10 +320,11 @@ int quiet(
double final_time,
double nu,
double delta,
double max_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,
@ -363,7 +366,7 @@ int quiet(
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);
} else if (algorithm==ALGORITHM_RKDP54) {
// only compute k1 on the first step
ns_step_rkdp54(u, adaptive_tolerance, adaptive_factor, max_delta, 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, 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);
@ -859,6 +862,7 @@ int ns_step_rkdp54(
double tolerance,
double factor,
double max_delta,
unsigned int adaptive_norm,
int K1,
int K2,
int N1,
@ -885,7 +889,6 @@ int ns_step_rkdp54(
){
int kx,ky;
double err,relative;
double sumu, sumU;
// 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)
@ -945,19 +948,54 @@ int ns_step_rkdp54(
// compute error
err=0;
sumu=0;
sumU=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_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)]);
// norm: |u_k|
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)]);
}
}
}
err=sqrt(err);
relative=pow((sqrt(sumu)+sqrt(sumU))/sumu, 1./3);
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++){
// 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)]);
}
}
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(err<relative*tolerance){
@ -979,7 +1017,7 @@ int ns_step_rkdp54(
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,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,false);
}
return 0;

View File

@ -33,13 +33,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, _Complex double* u0, _Complex double* g, bool irreversible, 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_norm, _Complex double* u0, _Complex double* g, bool irreversible, 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, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int algorithm, double print_freq, double starting_time, unsigned int nthreads, FILE* savefile, char* cmd_string, char* params_string, char* savefile_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_norm, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int algorithm, double print_freq, double starting_time, unsigned int nthreads, FILE* savefile, char* cmd_string, char* params_string, char* savefile_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, double starting_time, _Complex double* u0, _Complex double* g, bool irreversible, 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_norm, double starting_time, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int algorithm, unsigned int nthreads, FILE* savefile);
// initialize vectors for computation
@ -58,7 +58,7 @@ int ns_step_rk2( _Complex double* u, int K1, int K2, int N1, int N2, double nu,
// 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);
// Runge-Kutta-Dromand-Prince
int ns_step_rkdp54( _Complex double* u, double tolerance, double factor, double max_delta, 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_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);