From d81ad18618ba6b82727a65fc89c592b6bcee1ba5 Mon Sep 17 00:00:00 2001 From: Ian Jauslin Date: Mon, 18 Nov 2024 17:09:17 -0500 Subject: [PATCH] Rewrite cost function for adaptive step --- README.md | 8 +- docs/nstrophy_doc.tex | 85 ++++++------ src/constants.cpp | 9 +- src/main.c | 48 ++++--- src/navier-stokes.c | 307 +++++++++++++++++------------------------- src/navier-stokes.h | 18 ++- 6 files changed, 211 insertions(+), 264 deletions(-) diff --git a/README.md b/README.md index 810e3ef..1b87cd8 100644 --- a/README.md +++ b/README.md @@ -122,10 +122,10 @@ 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. +* `adaptive_cost`: cost function 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, `alpha` for alpha. * `keep_en_cst` (0 or 1, default 0): impose that the enstrophy is constant at each step (only really useful for the reversible equation). diff --git a/docs/nstrophy_doc.tex b/docs/nstrophy_doc.tex index 8ac0b90..f521e0e 100644 --- a/docs/nstrophy_doc.tex +++ b/docs/nstrophy_doc.tex @@ -115,18 +115,18 @@ In addition, several variable step methods are implemented: \end{itemize} In these adaptive step methods, two steps are computed at different orders: $\hat u_k^{(n)}$ and $\hat U_k^{(n)}$, the step size is adjusted at every step in such a way that the error is small enough: \begin{equation} - \|\hat u^{(n)}-\hat U^{(n)}\| + D(\hat u^{(n)},\hat U^{(n)}) <\epsilon_{\mathrm{target}} \end{equation} -for some given $\epsilon_{\mathrm{target}}$, set using the {\tt adaptive\_tolerance} parameter. -The choice of the norm matters, and will be discussed below. +for some given {\it cost function} $D$, and $\epsilon_{\mathrm{target}}$, set using the {\tt adaptive\_tolerance} parameter. +The choice of the cost function matters, and will be discussed below. If the error is larger than the target, then the step size is decreased. How this is done depends on the order of algorithm. -If the order is $q$ (here we mean the smaller of the two orders, so 4 for {\tt RKDP54} and {\tt RKF45} and 2 for {\tt RKBS32}), then we expect +If the order is $q$ (here we mean the smaller of the two orders, so 4 for {\tt RKDP54} and {\tt RKF45} and 2 for {\tt RKBS32}), then we expect (as long as the cost function is such that $D(\hat u,\hat u+\varphi)\sim\|\varphi\|$ in some norm) \begin{equation} - \|\hat u^{(n)}-\hat U^{(n)}\|=\delta_n^qC_n - . + D(\hat u^{(n)},\hat U^{(n)})=\delta_n^qC_n \end{equation} +for some number $C_n$. We wish to set $\delta_{n+1}$ so that \begin{equation} \delta_{n+1}^qC_n=\epsilon_{\mathrm{target}} @@ -135,7 +135,7 @@ so \begin{equation} \delta_{n+1} =\left(\frac{\epsilon_{\mathrm{target}}}{C_n}\right)^{\frac1q} - =\delta_n\left(\frac{\epsilon_{\mathrm{target}}}{\|\hat u^{(n)}-\hat U^{(n)}\|}\right)^{\frac1q} + =\delta_n\left(\frac{\epsilon_{\mathrm{target}}}{D(\hat u^{(n)},\hat U^{(n)})}\right)^{\frac1q} . \label{adaptive_delta} \end{equation} @@ -145,59 +145,49 @@ To be safe, we also set a maximal value for $\delta$ via the {\tt max\_delta} pa \bigskip \indent -The choice of the norm $\|\cdot\|$ matters. -It can be made by specifying the parameter {\tt adaptive\_norm}. +The choice of the cost function $D$ matters. +It can be made by specifying the parameter {\tt adaptive\_cost}. \begin{itemize} - \item A naive choice is to take $\|\cdot\|$ to be the normalized $L_1$ norm: + \item + For computations where the main focus is the enstrophy\-~(\ref{enstrophy}), one may want to set the cost function to the relative difference of the enstrophies: \begin{equation} - \|f\|:= - \frac1{\mathcal N}\sum_k|f_k| - ,\quad - \mathcal N:=\sum_k|\hat u_k^{(n)}-\hat u_k^{(n-1)}| + D(\hat u,\hat U):=\frac{|\mathcal En(\hat u)-\mathcal En(\hat U)|}{\mathcal En(\hat u)} . \end{equation} - This norm is selected by choosing {\tt adaptive\_norm=L1}. + This cost function is selected by choosing {\tt adaptive\_cost=enstrophy}. - \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 + \item + For computations where the main focus is the value of $\alpha$\-~(\ref{alpha}), one may want to set the cost function to the relative difference of $\alpha$: \begin{equation} - \|f\|:=\frac1{\mathcal N}\sum_k|f_k|k^{-3} + D(\hat u,\hat U):=\frac{|\alpha(\hat u)-\alpha(\hat U)|}{|\alpha(\hat u)|} + . + \end{equation} + This cost function is selected by choosing {\tt adaptive\_cost=alpha}. + + \item Alternatively, one my take $D$ to be the normalized $L_1$ norm: + \begin{equation} + D(\hat u,\hat U):= + \frac1{\mathcal N}\sum_k|\hat u_k-\hat U_k| ,\quad - \mathcal N:=\sum_k|\hat u_k^{(n)}-\hat u_k^{(n-1)}|k^{-3} + \mathcal N:=\sum_k|\hat u_k| + . + \end{equation} + This function is selected by choosing {\tt adaptive\_cost=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 cost function of the form + \begin{equation} + D(\hat u,\hat U):=\frac1{\mathcal N}\sum_k|\hat u_k-\hat U_k|k^{-3} + ,\quad + \mathcal N:=\sum_k|\hat u_k|k^{-3} \end{equation} or \begin{equation} - \|f\|:=\frac1{\mathcal N}\sum_k|f_k|k^{-\frac32} + D(\hat u,\hat U):=\frac1{\mathcal N}\sum_k|\hat u_k-\hat U_k|k^{-\frac32} ,\quad - \mathcal N:=\sum_k|\hat u_k^{(n)}-\hat u_k^{(n-1)}|k^{-\frac32} + \mathcal N:=\sum_k|\hat u_k|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}): - \begin{equation} - \|f\|:=\frac1{\mathcal N}\sqrt{\sum_k k^2|f_k|^2} - ,\quad - \mathcal N:=\frac{\sqrt{\sum_k k^2|\hat u_k^{(n)}|^2}+\sqrt{\sum_k k^2|\hat U_k^{(n)}|^2}}{\sum_k k^2|\hat u_k^{(n)}|^2} - . - \end{equation} - Doing so controls the error of the enstrophy through - \begin{equation} - \frac1{\mathcal N^2}|\mathcal En(\hat u)-\mathcal En(\hat U)|\equiv|\|\hat u\|^2-\|\hat U\|^2|\leqslant \|\hat u-\hat U\|(\|\hat u\|+\|\hat U\|) - \end{equation} - so - \begin{equation} - \frac1{\mathcal N^2} - |\mathcal En(\hat u)-\mathcal En(\hat U)|\leqslant - \|\hat u-\hat U\|\frac1{\mathcal N}\left(\sqrt{\sum_k k^2|\hat u_k|^2}+\sqrt{\sum_k k^2|\hat U_k|^2}\right) - \end{equation} - and thus - \begin{equation} - \frac{|\mathcal En(\hat u)-\mathcal En(\hat U)|}{\mathcal En(\hat u)}\leqslant - \|\hat u-\hat U\| - . - \end{equation} - This norm is selected by choosing {\tt adaptive\_norm=enstrophy}. + These cost functions are selected by choosing {\tt adaptive\_cost=k3} and {\tt adaptive\_cost=k32} respectively. \end{itemize} \subsubsection{Reality}. @@ -511,6 +501,7 @@ and so \alpha(\hat u) =\frac{\frac{L^2}{4\pi^2}\sum_k k^2\hat u_k^*\hat g_k+\sum_k|k|\hat u_k^*T(\hat u,k)}{\sum_kk^4|\hat u_k|^2} . + \label{alpha} \end{equation} Note that, by\-~(\ref{realu})-(\ref{realT}), \begin{equation} diff --git a/src/constants.cpp b/src/constants.cpp index 1cf1fd6..ca27947 100644 --- a/src/constants.cpp +++ b/src/constants.cpp @@ -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 diff --git a/src/main.c b/src/main.c index eebdaca..d1b045f 100644 --- a/src/main.c +++ b/src/main.c @@ -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(¶m_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(¶m_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); } } diff --git a/src/navier-stokes.c b/src/navier-stokes.c index 476ec4e..873ec59 100644 --- a/src/navier-stokes.c +++ b/src/navier-stokes.c @@ -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 || time0 ? -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(err0 ? -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(err0 ? -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(cost0 ? -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(err0 ? -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; diff --git a/src/navier-stokes.h b/src/navier-stokes.h index 9d9f989..79c0752 100644 --- a/src/navier-stokes.h +++ b/src/navier-stokes.h @@ -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);