From 251426aaaa81b76af12127b016f2b9b0594c70a9 Mon Sep 17 00:00:00 2001 From: Ian Jauslin Date: Tue, 13 Jun 2023 23:56:35 -0400 Subject: [PATCH] Select norm in parameter --- README.md | 5 +++ docs/nstrophy_doc.tex | 4 +++ src/constants.cpp | 4 +++ src/main.c | 43 ++++++++++++++++++++++++-- src/navier-stokes.c | 72 +++++++++++++++++++++++++++++++++---------- src/navier-stokes.h | 8 ++--- 6 files changed, 112 insertions(+), 24 deletions(-) diff --git a/README.md b/README.md index 8554692..c7eae36 100644 --- a/README.md +++ b/README.md @@ -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 diff --git a/docs/nstrophy_doc.tex b/docs/nstrophy_doc.tex index 24a5592..86f80da 100644 --- a/docs/nstrophy_doc.tex +++ b/docs/nstrophy_doc.tex @@ -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 diff --git a/src/constants.cpp b/src/constants.cpp index da593b8..53ec017 100644 --- a/src/constants.cpp +++ b/src/constants.cpp @@ -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 diff --git a/src/main.c b/src/main.c index 85b9c6d..34607b2 100644 --- a/src/main.c +++ b/src/main.c @@ -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){ diff --git a/src/navier-stokes.c b/src/navier-stokes.c index d04fe3b..ce699f6 100644 --- a/src/navier-stokes.c +++ b/src/navier-stokes.c @@ -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