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

delta_max and fix cabs2

Browse Source
This commit is contained in:
Ian Jauslin
2023-06-13 18:45:19 -04:00
parent 4ffbe1e978
commit 34b7a0c277
7 changed files with 166 additions and 20 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
README.md
View File

@@ -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

83
docs/nstrophy_doc.tex
View File

@@ -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

22
src/complex_tools.c Normal file
View File

@@ -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);
}

23
src/complex_tools.h Normal file
View File

@@ -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

21
src/main.c
View File

@@ -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){

26
src/navier-stokes.c
View File

@@ -17,13 +17,12 @@ limitations under the License.
#include "constants.cpp"
#include "io.h"
#include "navier-stokes.h"
#include "complex_tools.h"
#include <math.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
#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<relative*tolerance){
@@ -959,8 +965,8 @@ int ns_step_rkdp54(
u[klookup_sym(kx,ky,K2)]=tmp[klookup_sym(kx,ky,K2)];
}
}
// next delta to use in future steps
*next_delta=(*delta)*pow(relative*tolerance/err,1./5);
// next delta to use in future steps (do not exceed max_delta)
*next_delta=fmin(max_delta, (*delta)*pow(relative*tolerance/err,1./5));
// k1 in the next step will be this k4 (first same as last)
tmp=*k1;
@@ -971,7 +977,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,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,K1,K2,N1,N2,nu,delta,next_delta,L,g,fft1,fft2,ifft,k1,k2,k3,k4,k5,k6,tmp,irreversible,false);
}
return 0;

8
src/navier-stokes.h
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, _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, _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, _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, _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 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, 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, 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, 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);