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

Compare commits

base: software:7a0eb20452697ff65bc05efe5693c39f78832999
Branches Tags
software:master
...
compare: software:f68c179038f33bb7a762e0886c212a821695960b
Branches Tags
software:master

2 Commits
7a0eb20452 ... f68c179038

Author SHA1 Message Date
Ian Jauslin f68c179038 Error in enstrophy norm for rkf45 and rkbs32 2023-09-19 11:53:34 -04:00
Ian Jauslin d986f86243 Fix: no power 1/3 in enstrophy norm 2023-09-19 11:51:55 -04:00
2 changed files with 8 additions and 8 deletions
Show Stats Expand all files Collapse all files

6
docs/nstrophy_doc.tex
View File

@ -177,16 +177,16 @@ It can be made by specifying the parameter {\tt adaptive\_norm}.
\begin{equation}
\|f\|:=\frac1{\mathcal N}\sqrt{\sum_k k^2|f_k|^2}
,\quad
\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}
\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}
\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\|)
\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}
\mathcal N^2
\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}

10
src/navier-stokes.c
View File

@ -815,12 +815,12 @@ int ns_step_rkf45(
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((*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)]+(*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=pow((sqrt(sumu)+sqrt(sumU))/sumu, 1./3);
relative=(sqrt(sumu)+sqrt(sumU))/sumu;
}
else{
fprintf(stderr,"bug: unknown norm: %u, contact ian.jauslin@rutgers.edu\n", adaptive_norm);
@ -961,12 +961,12 @@ int ns_step_rkbs32(
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((*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(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=pow((sqrt(sumu)+sqrt(sumU))/sumu, 1./3);
relative=(sqrt(sumu)+sqrt(sumU))/sumu;
}
else{
fprintf(stderr,"bug: unknown norm: %u, contact ian.jauslin@rutgers,edu\n", adaptive_norm);
@ -1143,7 +1143,7 @@ int ns_step_rkdp54(
}
}
err=sqrt(err);
relative=pow((sqrt(sumu)+sqrt(sumU))/sumu, 1./3);
relative=(sqrt(sumu)+sqrt(sumU))/sumu;
}
else{
fprintf(stderr,"bug: unknown norm: %u, contact ian.jauslin@rutgers,edu\n", adaptive_norm);