From 34b7a0c277f2dfeec07197650a354030d5763ec8 Mon Sep 17 00:00:00 2001 From: Ian Jauslin Date: Tue, 13 Jun 2023 18:45:19 -0400 Subject: [PATCH] delta_max and fix cabs2 --- README.md | 3 ++ docs/nstrophy_doc.tex | 83 +++++++++++++++++++++++++++++++++++++++++++ src/complex_tools.c | 22 ++++++++++++ src/complex_tools.h | 23 ++++++++++++ src/main.c | 21 +++++++---- src/navier-stokes.c | 26 ++++++++------ src/navier-stokes.h | 8 ++--- 7 files changed, 166 insertions(+), 20 deletions(-) create mode 100644 src/complex_tools.c create mode 100644 src/complex_tools.h diff --git a/README.md b/README.md index b605b22..8554692 100644 --- a/README.md +++ b/README.md @@ -105,6 +105,9 @@ should be a `;` sperated list of `key=value` pairs. The possible keys are * `adaptive_factor` (double, default 0.9): when using the RKF45 method, the step gets adjusted by `factor*delta*(tolerance/error)^(1/5)`. +* `max_delta` (double, default 1e-2): when using the adaptive step, do not + exceet `delta_max`. + # Interrupting/resuming the computation diff --git a/docs/nstrophy_doc.tex b/docs/nstrophy_doc.tex index d05d23f..ba15feb 100644 --- a/docs/nstrophy_doc.tex +++ b/docs/nstrophy_doc.tex @@ -103,6 +103,88 @@ We truncate the Fourier modes and assume that $\hat u_k=0$ if $|k_1|>K_1$ or $|k \end{equation} \bigskip +\point{\bf Runge-Kutta methods}. +To solve the equation numerically, we will use Runge-Kutta methods, which compute an approximate value $\hat u_k^{(n)}$ for $\hat u_k(t_n)$. +{\tt nstrophy} supports the 4th order Runge-Kutta ({\tt RK4}) and 2nd order Runge-Kutta ({\tt RK2}) algorithms. +In addition, several variable step methods are implemented: +\begin{itemize} + \item the Runge-Kutta-Dormand-Prince method ({\tt RKDP54}), which is of 5th order, and adjusts the step by comparing to a 4th order method; + \item the Runge-Kutta-Fehlberg method ({\tt RKF45}), which is of 4th order, and adjusts the step by comparing to a 5th order method; + \item the Runge-Kutta-Bogacki-Shampine method ({\tt RKBS32}), which is of 3d order, and adjusts the step by comparing to a 2nd order method. +\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)}\| + <\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. +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 +\begin{equation} + \|\hat u^{(n)}-\hat U^{(n)}\|=\delta_n^qC_n + . +\end{equation} +We wish to set $\delta_{n+1}$ so that +\begin{equation} + \delta_{n+1}^qC_n=\epsilon_{\mathrm{target}} +\end{equation} +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} + . + \label{adaptive_delta} +\end{equation} +(Actually, to be safe and ensure that $\delta$ decreases sufficiently, we multiply this by a safety factor that can be set using the {\tt adaptive\_factor} parameter.) +If the error is smaller than the target, we increase $\delta$ using\-~(\ref{adaptive_delta}) (without the safety factor). +To be safe, we also set a maximal value for $\delta$ via the {\tt max\_delta} parameter. +\bigskip + +\indent +The choice of the norm $\|\cdot\|$ matters. +\begin{itemize} + \item A naive choice is to take $\|\cdot\|$ to be the normalized $L_1$ norm: + \begin{equation} + \|f\|:= + \frac1{\mathcal N}\sum_k|f_k| + ,\quad + \mathcal N:=\sum_k|\hat u_k^{(n)}-\hat u_k^{(n-1)}| + . + \end{equation} + + \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} + \|f\|:=\frac1{\mathcal N}\sum_k|f_k|k^{-3} + ,\quad + \mathcal N:=\sum_k|\hat u_k^{(n)}-\hat u_k^{(n-1)}|k^{-3} + \end{equation} + or + \begin{equation} + \|f\|:=\frac1{\mathcal N}\sum_k|f_k|k^{-\frac32} + ,\quad + \mathcal N:=\sum_k|\hat u_k^{(n)}-\hat u_k^{(n-1)}|k^{-\frac32} + \end{equation} + are sensible choices. + + \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:=\sqrt{\sum_k k^2|\hat u_k^{(n)}|^2} + . + \end{equation} + Doing so controls the error of the enstrophy through + \begin{equation} + \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} +\end{itemize} +\bigskip + \point{\bf Reality}. Since $U$ is real, $\hat U_{-k}=\hat U_k^*$, and so \begin{equation} @@ -271,6 +353,7 @@ and +e^{-\frac{4\pi^2}{L^2}\nu t}\sqrt{E(0)} \right)^2 . + \label{enstrophy} \end{equation} \bigskip diff --git a/src/complex_tools.c b/src/complex_tools.c new file mode 100644 index 0000000..1afbc62 --- /dev/null +++ b/src/complex_tools.c @@ -0,0 +1,22 @@ +/* +Copyright 2017-2023 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +#include "complex_tools.h" + +// magnitude squared +double cabs2(_Complex double x){ + return (__real__ x)*(__real__ x)+(__imag__ x)*(__imag__ x); +} diff --git a/src/complex_tools.h b/src/complex_tools.h new file mode 100644 index 0000000..c9802d5 --- /dev/null +++ b/src/complex_tools.h @@ -0,0 +1,23 @@ +/* +Copyright 2017-2023 Ian Jauslin + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +#ifndef COMPLEXTOOLS_H +#define COMPLEXTOOLS_H + +// magnitude squared +double cabs2(_Complex double x); + +#endif diff --git a/src/main.c b/src/main.c index 011dcae..85b9c6d 100644 --- a/src/main.c +++ b/src/main.c @@ -47,6 +47,7 @@ typedef struct nstrophy_parameters { double L; double adaptive_tolerance; double adaptive_factor; + double max_delta; double print_freq; int seed; double starting_time; @@ -170,16 +171,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, 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, 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, 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, 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.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.starting_time, u0, g, parameters.irreversible, parameters.algorithm, nthreads, savefile); } else if(command==0){ fprintf(stderr, "error: no command specified\n"); @@ -264,13 +265,13 @@ int print_params( fprintf(file,", algorithm=RK2"); break; case ALGORITHM_RKF45: - fprintf(file,", algorithm=RKF45, tolerance=%.15e, factor=%.15e",parameters.adaptive_tolerance, parameters.adaptive_factor); + fprintf(file,", algorithm=RKF45, tolerance=%.15e",parameters.adaptive_tolerance); break; case ALGORITHM_RKDP54: - fprintf(file,", algorithm=RKDP54, tolerance=%.15e, factor=%.15e",parameters.adaptive_tolerance, parameters.adaptive_factor); + fprintf(file,", algorithm=RKDP54, tolerance=%.15e",parameters.adaptive_tolerance); break; case ALGORITHM_RKBS32: - fprintf(file,", algorithm=RKBS32, tolerance=%.15e, factor=%.15e",parameters.adaptive_tolerance, parameters.adaptive_factor); + fprintf(file,", algorithm=RKBS32, tolerance=%.15e",parameters.adaptive_tolerance); break; default: fprintf(file,", algorithm=RK4"); @@ -397,6 +398,7 @@ int read_params( parameters->L=2*M_PI; parameters->adaptive_tolerance=1e-11; parameters->adaptive_factor=0.9; + parameters->max_delta=1e-2; parameters->final_time=100000; parameters->print_freq=1; parameters->starting_time=0; @@ -597,6 +599,13 @@ int set_parameter( return(-1); } } + else if (strcmp(lhs,"max_delta")==0){ + ret=sscanf(rhs,"%lf",&(parameters->max_delta)); + if(ret!=1){ + fprintf(stderr, "error: parameter 'max_delta' should be a double\n got '%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 186864e..996b390 100644 --- a/src/navier-stokes.c +++ b/src/navier-stokes.c @@ -17,13 +17,12 @@ limitations under the License. #include "constants.cpp" #include "io.h" #include "navier-stokes.h" +#include "complex_tools.h" #include #include #include #include -#define CABS2(x) ((__real__ x) * (__real__ x) + (__imag__ x) * (__imag__ x)) - // compute solution as a function of time int uk( int K1, @@ -36,6 +35,7 @@ int uk( double L, double adaptive_tolerance, double adaptive_factor, + double max_delta, _Complex double* u0, _Complex double* g, bool irreversible, @@ -94,7 +94,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, 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, 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); @@ -144,6 +144,7 @@ int enstrophy( double L, double adaptive_tolerance, double adaptive_factor, + double max_delta, _Complex double* u0, _Complex double* g, bool irreversible, @@ -203,7 +204,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, 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, 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); @@ -317,6 +318,7 @@ int quiet( double final_time, double nu, double delta, + double max_delta, double L, double adaptive_tolerance, double adaptive_factor, @@ -361,7 +363,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, 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, 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); @@ -856,6 +858,7 @@ int ns_step_rkdp54( _Complex double* u, double tolerance, double factor, + double max_delta, int K1, int K2, int N1, @@ -946,10 +949,13 @@ int ns_step_rkdp54( 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)])); - relative+=(kx*kx+ky*ky)*(CABS2(tmp[klookup_sym(kx,ky,K2)]-u[klookup_sym(kx,ky,K2)])); + 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)])); + //relative+=(kx*kx+ky*ky)*(CABS2(tmp[klookup_sym(kx,ky,K2)]-u[klookup_sym(kx,ky,K2)])); + relative+=(kx*kx+ky*ky)*(cabs2(u[klookup_sym(kx,ky,K2)])); } } + err=sqrt(err); + relative=sqrt(relative); // compare relative error with tolerance if(err