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

Proper normalization for enstrophy norm

Browse Source
This commit is contained in:
Ian Jauslin 2023-06-13 19:01:37 -04:00
parent 34b7a0c277
commit 3b58896e3b
2 changed files with 18 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
docs/nstrophy_doc.tex
View File

@ -174,12 +174,23 @@ The choice of the norm $\|\cdot\|$ matters.
\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}
\mathcal N:=\left(\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}\right)^{\frac13}
.
\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}
so
\begin{equation}
\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}
\end{itemize}

10
src/navier-stokes.c
View File

@ -885,6 +885,7 @@ 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)
@ -944,18 +945,19 @@ int ns_step_rkdp54(
// compute error
err=0;
relative=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)]));
//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)]));
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(relative);
relative=pow((sqrt(sumu)+sqrt(sumU))/sumu, 1./3);
// compare relative error with tolerance
if(err<relative*tolerance){