sandbox/ghigo/artery1D/thacker.c
Thacker solution
We solve, using the 1D blood flow equations, the inviscid pressure oscillations in a parabolic aneurysm. This test case is greatly inspired by the solutions of Thacker and Sampson, well-known in the shallow water community.
Analytic solution
We choose the following spatial variations of the neutral cross-sectional area A_0 and the arterial wall rigidity K^\prime in the domain x\in \left[-a,\, a \right]:
\displaystyle \left\{ \begin{aligned} & R_0 = \bar{R_0}\left[1 + \Delta \mathcal{G} \left[1 - \frac{x^2}{a^2} \right] \right] && \text{ with } \bar{R_0} >0 ,\, a >0 ,\, \Delta \mathcal{G} > -1 \\ & K^\prime = \mathrm{cst} && \text{ with } K^\prime >0, \end{aligned} \right.
where \Delta \mathcal{G} > 0. Simple analysis allows us to obtain the following analytic solution:
\displaystyle \left\{ \begin{aligned} & U_0 = U \sin \left( \frac{t}{\tau} \right) && \text{ with } U = \mathrm{cst}\\ & p = - \frac{1}{4} \rho U^2 \cos \left( 2 \frac{t}{\tau} \right) - \rho \frac{x}{\tau} U \cos \left( \frac{t}{\tau} \right), \end{aligned} \right.
where:
\displaystyle \tau = \lvert \Delta \mathcal{G} \omega^2 \rvert^{-\frac{1}{2}} .
Finally, we choose U with respect to the following nonlinear stability arguments, namely that the radius R remains positive and that the flow remains subcritical.
Code
#include "grid/cartesian1D.h"
#include "bloodflow.h"
double R0 = 1.;
double K0 = 1.e4;
double BTH = 4.;
double DELTATH = 0.1;
double ATH = 1.;
double KTH = 0.;
double CTH = 0.;
double WTH = 0.;
double TAUTH = 0.;
double TTH = 0.;
scalar eq[];
double eq1, eq2, eqmax;
int ne;
int ihr = 0;
double celerity (double a, double k) {
return sqrt(0.5*k*sqrt(a));
}
double r0thacker(double x, double r0, double b, double delta) {
return r0*(1. + delta*(1. - pow(x/b, 2.)));
}
double pthacker(double x, double t, double r0, double rho,
double tau, double ath) {
return -0.25*rho*pow(ath, 2.)*cos(2.*t/tau) - rho*ath*x/tau*cos(t/tau);
}
double rthacker(double x, double t, double r0, double b,
double delta, double k, double rho, double tau, double ath) {
return 1./k*pthacker(x, t, r0, rho, tau, ath) + r0thacker(x, r0, b, delta);
}
double athacker(double x, double t, double r0, double b,
double delta, double k, double rho, double tau, double ath) {
return pi*pow(rthacker(x, t, r0, b, delta, k, rho, tau, ath), 2.);
}
double uthacker(double t, double tau, double ath) {
return ath*sin(t/tau);
}
double qthacker(double x, double t, double r0, double b,
double delta, double k, double rho, double tau, double ath) {
return uthacker(t, tau, ath)*athacker(x, t, r0, b, delta, k, rho, tau, ath);
}
int main() {
origin (-BTH, 0.);
L0 = 2.*BTH;
KTH = K0*sqrt(pi);
CTH = sqrt(0.5*KTH*R0);
WTH = 2.*CTH/BTH;
TAUTH = 1./sqrt(abs(DELTATH*pow(WTH, 2.)));
TTH = 2.*pi/sqrt(abs(DELTATH*pow(WTH, 2.)));
ihr = 0;
for (ihr = 0; ihr <= 2; ihr ++) {
for (N = 32; N <= 256; N *= 2) {
eq1 = eq2 = eqmax = 0.;
ne = 0;
run();
printf ("eq%d, %d, %g, %g, %g\n", ihr, N, eq1/ne,
sqrt(eq2/ne), eqmax);
}
}
return 0;
}
q[left] = dirichlet (qthacker (-BTH, t, R0, BTH, DELTATH, KTH, 1., TAUTH, ATH));
q[right] = dirichlet (qthacker (BTH, t, R0, BTH, DELTATH, KTH, 1., TAUTH, ATH));
event defaults (i=0) {
gradient = order1;
if (ihr == 0) riemann = hll_hr;
else if (ihr == 1) riemann = hll_hrls;
else if (ihr == 2) riemann = hll_glu;
}
event init (i=0) {
foreach() {
k[] = K0;
a0[] = pi*pow(r0thacker (x, R0, BTH, DELTATH), 2.);
a[] = athacker (x, 0., R0, BTH, DELTATH, KTH, 1., TAUTH, ATH);
q[] = qthacker (x, 0., R0, BTH, DELTATH, KTH, 1., TAUTH, ATH);
}
}
event field (t = {0., 0.1*0.422653, 0.3*0.422653, 0.6*0.422653, 0.8*0.422653}) {
if (N == 128) {
foreach() {
fprintf (stderr, "%d, %g, %.6f, %.6f, %.6f, %.6f\n", ihr, x,
a0[]/(pi*R0*R0), k[]/K0, a[]-a0[], q[]);
}
fprintf (stderr, "\n\n");
}
}
event error (i++) {
ne++;
foreach() {
eq[] = q[] - qthacker (x, t, R0, BTH, DELTATH,
KTH, 1., TAUTH, ATH);
}
norm nq = normf (eq);
eq1 += nq.avg;
eq2 += nq.rms*nq.rms;
if (nq.max > eqmax)
eqmax = nq.max;
}
event end (t = 1.) {
printf ("#Thacker test case completed\n");
}
Plots
reset
mycolors = "dark-blue red sea-green dark-violet orange brown4"
mypoints = '1 2 3 4 6'
set datafile separator ','
set style line 1 lt -1 lw 2 lc 'black'
set output 'A0K.png'
set xlabel 'x [cm]'
set ylabel 'Dimensionless arterial properties'
set key top right
plot '< grep -e "^0" -e "^$" log' i 0 u 2:3 w l ls 1 lc rgb word(mycolors, 1) t 'A_0/A_0(0)', \
'< grep -e "^0" -e "^$" log' i 0 u 2:4 w l ls 1 lc rgb word(mycolors, 2) t 'K/K(0)'
r0 = 1.
k0 = 1.e4
bth = 4.
deltath = 0.1
ath = 1.
kth = k0*sqrt(pi)
cth = sqrt(0.5*kth*r0)
wth = 2.*cth/bth
tauth = 1./sqrt(abs(deltath*(wth)**2.))
tth = 2.*pi/sqrt(abs(deltath*(wth)**2.))
r0thacker(x) = r0*(1. + deltath*(1. - (x/bth)**2.))
a0thacker(x) = pi*r0thacker(x)**2.
pthacker(x, t) = -0.25*(ath)**2.*cos(2.*t/tauth) - ath*x/tauth*cos(t/tauth)
rthacker(x, t) = 1./kth*pthacker(x, t) + r0thacker(x)
athacker(x, t) = pi*rthacker(x, t)**2.
uthacker(t) = ath*sin(t/tauth)
qthacker(x, t) = uthacker(t)*athacker(x, t)
set output 'AmA0.png'
set xlabel 'x [cm]'
set ylabel 'A-A_0 [cm^2]'
set key top center
plot athacker(x, 0) - a0thacker(x) w l ls 1 t 'Analytic', \
athacker(x, 0.1*0.422653) - a0thacker(x) w l ls 1 notitle, \
athacker(x, 0.3*0.422653) - a0thacker(x) w l ls 1 notitle, \
athacker(x, 0.6*0.422653) - a0thacker(x) w l ls 1 notitle, \
athacker(x, 0.8*0.422653) - a0thacker(x) w l ls 1 notitle, \
'< grep -e "^0" -e "^$" log' i 0 u 2:5 every 2 w p lt word(mypoints, 1) lc rgb word(mycolors, 1) ps 2 lw 1 t 't=0', \
'< grep -e "^0" -e "^$" log' i 1 u 2:5 every 2 w p lt word(mypoints, 2) lc rgb word(mycolors, 2) ps 2 lw 1 t 't=0.1T', \
'< grep -e "^0" -e "^$" log' i 2 u 2:5 every 2 w p lt word(mypoints, 3) lc rgb word(mycolors, 3) ps 2 lw 1 t 't=0.3T', \
'< grep -e "^0" -e "^$" log' i 3 u 2:5 every 2 w p lt word(mypoints, 4) lc rgb word(mycolors, 4) ps 2 lw 1 t 't=0.6T', \
'< grep -e "^0" -e "^$" log' i 4 u 2:5 every 2 w p lt word(mypoints, 5) lc rgb word(mycolors, 5) ps 2 lw 1 t 't=0.8T'
{0,, 0.1,, 0.3,, 0.6,, 0.8}
set output 'Q.png'
set xlabel 'x [cm]'
set ylabel 'Q [cm^3/s]'
set key center center
plot qthacker(x, 0) w l ls 1 t 'Analytic', \
qthacker(x, 0.1*0.422653) w l ls 1 notitle, \
qthacker(x, 0.3*0.422653) w l ls 1 notitle, \
qthacker(x, 0.6*0.422653) w l ls 1 notitle, \
qthacker(x, 0.8*0.422653) w l ls 1 notitle, \
'< grep -e "^0" -e "^$" log' i 0 u 2:6 every 2 w p lt word(mypoints, 1) lc rgb word(mycolors, 1) ps 2 lw 1 t 't=0', \
'< grep -e "^0" -e "^$" log' i 1 u 2:6 every 2 w p lt word(mypoints, 2) lc rgb word(mycolors, 2) ps 2 lw 1 t 't=0.1T', \
'< grep -e "^0" -e "^$" log' i 2 u 2:6 every 2 w p lt word(mypoints, 3) lc rgb word(mycolors, 3) ps 2 lw 1 t 't=0.3T', \
'< grep -e "^0" -e "^$" log' i 3 u 2:6 every 2 w p lt word(mypoints, 4) lc rgb word(mycolors, 4) ps 2 lw 1 t 't=0.6T', \
'< grep -e "^0" -e "^$" log' i 4 u 2:6 every 2 w p lt word(mypoints, 5) lc rgb word(mycolors, 5) ps 2 lw 1 t 't=0.8T'
{0,, 0.1,, 0.3,, 0.6,, 0.8}
set xtics 32,2,256
set logscale
ftitle(a,b) = sprintf('Order %4.2f', -b)
f1(x) = a1 + b1*x
f2(x) = a2 + b2*x
f3(x) = a3 + b3*x
fit f1(x) '< grep ^eq0 out' u (log($2)):(log($4)) via a1, b1
fit f2(x) '< grep ^eq1 out' u (log($2)):(log($4)) via a2, b2
fit f3(x) '< grep ^eq2 out' u (log($2)):(log($4)) via a3, b3
set output 'errQL2.png'
set xlabel 'N'
set ylabel '|Q|_2'
set key top right
set cbrange[1e-10:1e10]
plot '< grep ^eq0 out' u 2:4 w p lt word(mypoints, 1) lc rgb word(mycolors, 1) ps 3 lw 2 t 'HR, '.ftitle(a1, b1), \
exp (f1(log(x))) ls 1 notitle, \
'< grep ^eq1 out' u 2:4 w p lt word(mypoints, 2) lc rgb word(mycolors, 2) ps 3 lw 2 t 'HR LS, '.ftitle(a2, b2), \
exp (f2(log(x))) ls 1 notitle, \
'< grep ^eq2 out' u 2:4 w p lt word(mypoints, 3) lc rgb word(mycolors, 3) ps 3 lw 2 t 'GLU, '.ftitle(a3, b3), \
exp (f3(log(x))) ls 1 notitle
fit f1(x) '< grep ^eq0 out' u (log($2)):(log($3)) via a1, b1
fit f2(x) '< grep ^eq0 out' u (log($2)):(log($4)) via a2, b2
fit f3(x) '< grep ^eq0 out' u (log($2)):(log($5)) via a3, b3
set output 'errQHR.png'
set xlabel 'N'
set ylabel 'Relative error norms for HR'
set key top right
plot '< grep ^eq0 out' u 2:3 w p lt word(mypoints, 1) lc rgb word(mycolors, 1) ps 3 lw 2 t '|Q|_1, '.ftitle(a1, b1), \
exp (f1(log(x))) ls 1 notitle, \
'< grep ^eq0 out' u 2:4 w p lt word(mypoints, 2) lc rgb word(mycolors, 2) ps 3 lw 2 t '|Q|_2, '.ftitle(a2, b2), \
exp (f2(log(x))) ls 1 notitle, \
'< grep ^eq0 out' u 2:5 w p lt word(mypoints, 3) lc rgb word(mycolors, 3) ps 3 lw 2 t '|Q|_{max}, '.ftitle(a3, b3), \
exp (f3(log(x))) ls 1 notitle