tourniquet-adapt.c

Inviscid tourniquet
Results for first order

Inviscid tourniquet

We solve here the inviscid 1D blood flow equations in a straight artery. The artery is initially inflated on its left half, mimicking the action of a tourniquet. At t>0, the vessel relaxes towards its steady-state at rest. The solution of the flow of blood generated by this elastic relaxation is obtained numerically and compared to the solution of the corresponding Riemann problem.

Analytic solution of the corresponding Riemann problem

The initial conditions are: \displaystyle A(t=0,\, x) = \left\{ \begin{aligned} & A_L \quad \mathrm{if} \: x < x_m \\ & A_R \quad \mathrm{if} \: x >= x_m \end{aligned} \right. , \qquad Q(t=0,\, x) = 0, with A_L > A_R and x_m = 0. Using the method of characteristic, we obtain the following analytic solution: \displaystyle \begin{aligned} & \mathrm{if} \: x < x_A = - c_L t \quad \left\{ \begin{aligned} & A(t,\, x) = A_L \\ & U(t,\, x) = 0 \end{aligned} \right. \\ & \mathrm{if} \: x_A <= x < x_B = (u_M - c_M)t \quad \left\{ \begin{aligned} & U(t,\, x) = \frac{4}{5}\frac{x}{t} + \frac{4}{5}c_L \\ & c(t,\, x) = -\frac{1}{5}\frac{x}{t} + \frac{4}{5}c_L \end{aligned} \right. \\ & \mathrm{if} \: x_B <= x < x_C = \frac{A_M U_M}{A_M - A_R} t \quad \left\{ \begin{aligned} & A(t,\, x) = A_M \\ & U(t,\, x) = U_M \end{aligned} \right. \\ & \mathrm{if} \: x_C <= x \quad \left\{ \begin{aligned} & A(t,\, x) = A_R \\ & U(t,\, x) = 0 \end{aligned} \right. . \end{aligned} The variables A_M and U_M are obtain by solving the following system iteratively: \displaystyle \left\{ \begin{aligned} & U_M + 4 c_M = 4 c_L\\ & Q_M = \frac{A_M U_M}{A_M - A_R} \left[A_M - A_R\right]\\ & \left[ \frac{Q_M^2}{A_M} + \frac{K}{3\rho} A_M^{\frac{3}{2}} \right] - \frac{K}{3 \rho} A_R^{\frac{3}{2}} = \frac{A_M U_M}{A_M - A_R} Q_M . \end{aligned} \right.

double celerity (double a, double k)
{
  return sqrt(0.5*k*sqrt(a));
}

double analytic_a (double t, double x, double dr,
		   double a0, double k0)
{
  double al = a0*pow(1 + dr, 2.);
  double cl = celerity(al, k0);

  double ar = a0;

  double am = 3.459578046858399;
  double cm = celerity(am, k0);
  double um = 9.192473939896399;
  double s = am*um/(am-ar);

  double xa = -cl*t;
  double xb = (um - cm)*t;
  double xc = s*t;

  if (x <= xa) return al;
  else if (xa < x && x <= xb) {
    return pow(2./k0, 2.)*pow( -1./5.*x/t + 4./5.*cl, 4.);
  }
  else if (xb < x && x <= xc) return am;
  else return ar;	     
}

double analytic_u (double t, double x, double dr,
		  double a0, double k0)
{
  double al = a0*pow(1 + dr, 2.);
  double cl = celerity(al, k0);

  double ar = a0;
  
  double am = 3.459578046858399;
  double cm = celerity(am, k0);
  double um = 9.192473939896399;
  double s = am*um/(am-ar);

  double xa = -cl*t;
  double xb = (um - cm)*t;
  double xc = s*t;

  if (x <= xa) return 0.;
  else if (xa < x && x <= xb) return  4./5.*x/t + 4./5.*cl;
  else if (xb < x && x <= xc) return um;
  else return 0.;	     
}

double analytic_q (double t, double x, double dr,
		  double a0, double k0)
{
  return analytic_a(t, x, dr, a0, k0)*analytic_u(t, x, dr, a0, k0);
}

We choose here to use an adaptive cartesian grid to capture the sharp profile of the perturbation of the artery.

#include "grid/bitree.h"
#include "../bloodflow-hr.h"

#define lmin (5) // N = 32
int lmax;
#define cmax (1.e-3)

We define the artery’s geometrical and mechanical properties.

#define R0 (1.)
#define XM (0.)
#define DR (1.e-1)
#define shape(x) ((x) < (XM) ? (1. + (DR)) : 1.)

#define K0 (1.e4)

int main() {

The domain is 10.

  L0 = 10.;
  size (L0);
  origin (-(L0)/2.);
  
  DT = 1.e-5;

We run the computation for maximum level of refinements lmax.

  for (lmax = 8; lmax <= 14; lmax++) {
    N = (1 << (lmax));
    init_grid (N);
    run();
  }
}

Boundary conditions

We impose homogeneous Neumann boundary conditions on all variables.

a[left] = neumann (0.);
q[left] = neumann (0.);

a[right] = neumann (0.);
q[right] = neumann (0.);

Defaults conditions

event defaults (i = 0)
{
  gradient = zero;
}

Initial conditions

event init (i = 0) {

We initialize the variables k, zb, a and q.

  foreach() {
    k[] = (K0);
    zb[] = k[]*sqrt (pi)*(R0);
    a[] = sq (zb[]/k[]*(shape (x)));
    q[] = 0.;
  }
}

Post-processing

We output the computed fields.

event field (t = {0., 0.01, 0.02, 0.03, 0.04}) {

  if (lmax == 14) {

    char name[80];
    sprintf (name, "fields-%.2f-pid-%d.dat", t, pid());
    FILE * ff = fopen (name, "w");
    
    foreach()
      fprintf (ff, "%g %g %g %g %g %g %g %d\n",
	       x, k[], sq (zb[]/k[]),
	       (analytic_a (t,x,(DR),(pi*sq ((R0))),(K0)))/(pi*sq ((R0))),
	       a[]/(pi*sq ((R0))),
	       (analytic_q (t,x,(DR),(pi*sq ((R0))),(K0))),
	       q[],
	       level
	       );
  }
}

Next, we compute the spatial error for the cross-sectional area a and the flow rate q.

event error (t = 0.04)
{
  scalar err_a[], err_q[];
  foreach() {
    err_a[] = fabs (a[] - analytic_a (t, x, (DR), (pi*sq ((R0))), (K0)));
    err_q[] = fabs (q[] - analytic_q (t, x, (DR), (pi*sq ((R0))), (K0)));
  }
  boundary ((scalar *) {err_a, err_q});
  
  norm na = normf (err_a);
  norm nq = normf (err_q);
  
  fprintf (ferr, "%d %g %g %g %g %g %g\n",
	   (lmax),
	   na.avg, na.rms, na.max,
	   nq.avg, nq.rms, nq.max);
}

Mesh adaptation

event adapt (i++)
{
  adapt_wavelet ({a,q}, (double[]){(cmax),(cmax)},
		 maxlevel = (lmax), minlevel = (lmin));
}

End of simulation

event stop_run (t = 0.05)
{
  return 0;
}

Results for first order

Mesh adaptation

We first the distribution of the mesh levels throughout the artery at t={0, 0.01, 0.02, 0.03, 0.04} for lmax=14.

reset
set xlabel 'x'
set ylabel 'l'
set yrange[:15]
plot '< cat fields-0.00-pid-*' u 1:8 w l lw 2 lc rgb "blue" t 't=0', \
     '< cat fields-0.01-pid-*' u 1:8 w l lw 2 lc rgb "red" t 't=0.01', \
     '< cat fields-0.02-pid-*' u 1:8 w l lw 2 lc rgb "sea-green" t 't=0.02', \
     '< cat fields-0.03-pid-*' u 1:8 w l lw 2 lc rgb "coral" t 't=0.03', \
     '< cat fields-0.04-pid-*' u 1:8 w l lw 2 lc rgb "dark-violet" t 't=0.04'

levels for lmax=14. (script)

Cross-sectional area and flow rate

We then plot the spatial evolution of the cross-sectional area a at t={0, 0.01, 0.02, 0.03, 0.04} for lmax=14.

reset
set xlabel 'x'
set ylabel 'a/a_0'
plot '< cat fields-0.00-pid-*' u 1:4 w l lw 3 lc rgb "black" t 'analytic', \
     '< cat fields-0.01-pid-*' u 1:4 w l lw 3 lc rgb "black" notitle, \
     '< cat fields-0.02-pid-*' u 1:4 w l lw 3 lc rgb "black" notitle, \
     '< cat fields-0.03-pid-*' u 1:4 w l lw 3 lc rgb "black" notitle, \
     '< cat fields-0.04-pid-*' u 1:4 w l lw 3 lc rgb "black" notitle, \
     '< cat fields-0.00-pid-*' u 1:5 w l lw 2 lc rgb "blue" t 't=0', \
     '< cat fields-0.01-pid-*' u 1:5 w l lw 2 lc rgb "red" t 't=0.01', \
     '< cat fields-0.02-pid-*' u 1:5 w l lw 2 lc rgb "sea-green" t 't=0.02', \
     '< cat fields-0.03-pid-*' u 1:5 w l lw 2 lc rgb "coral" t 't=0.03', \
     '< cat fields-0.04-pid-*' u 1:5 w l lw 2 lc rgb "dark-violet" t 't=0.04'

a/a_0 for lmax=14. (script)

reset
set key bottom right
set xlabel 'x'
set ylabel 'q'
plot '< cat fields-0.00-pid-*' u 1:6 w l lw 3 lc rgb "black" t 'analytic', \
     '< cat fields-0.01-pid-*' u 1:6 w l lw 3 lc rgb "black" notitle, \
     '< cat fields-0.02-pid-*' u 1:6 w l lw 3 lc rgb "black" notitle, \
     '< cat fields-0.03-pid-*' u 1:6 w l lw 3 lc rgb "black" notitle, \
     '< cat fields-0.04-pid-*' u 1:6 w l lw 3 lc rgb "black" notitle, \
     '< cat fields-0.00-pid-*' u 1:7 w l lw 2 lc rgb "blue" t 't=0', \
     '< cat fields-0.01-pid-*' u 1:7 w l lw 2 lc rgb "red" t 't=0.01', \
     '< cat fields-0.02-pid-*' u 1:7 w l lw 2 lc rgb "sea-green" t 't=0.02', \
     '< cat fields-0.03-pid-*' u 1:7 w l lw 2 lc rgb "coral" t 't=0.03', \
     '< cat fields-0.04-pid-*' u 1:7 w l lw 2 lc rgb "dark-violet" t 't=0.04'

q for lmax=14. (script)

Convergence

Finally, we plot the evolution of the error for the cross-sectional area a and the flow rate q with the number of cells N_{max}. We compare the results with those obtained with a uniform grid.

reset
set xlabel 'N_{max}'
set ylabel 'L_1(a),L_2(a),L_{max}(a)'
set format y '%.1e'
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) '../tourniquet/log' u (log($1)):(log($2)) via a1, b1
fit f2(x) '../tourniquet/log' u (log($1)):(log($3)) via a2, b2
fit f3(x) '../tourniquet/log' u (log($1)):(log($4)) via a3, b3

f11(x) = a11 + b11*x
f22(x) = a22 + b22*x
f33(x) = a33 + b33*x
fit f11(x) 'log' u (log(2**$1)):(log($2)) via a11, b11
fit f22(x) 'log' u (log(2**$1)):(log($3)) via a22, b22
fit f33(x) 'log' u (log(2**$1)):(log($4)) via a33, b33

f111(x) = a111 + b111*x
f222(x) = a222 + b222*x
fit [x=2^10:] f111(x) 'log' u (log(2**$1)):(log($2)) via a111, b111
fit [x=2^10:] f222(x) 'log' u (log(2**$1)):(log($3)) via a222, b222

plot '../tourniquet/log' u 1:2 w p pt 6 ps 1.5 lc rgb "blue" t '|a|_1, '.ftitle(a1, b1), \
     exp (f1(log(x))) ls 1 lc rgb "red" notitle, \
     'log' u (2**$1):2 w p pt 6 ps 1.5 lc rgb "sea-green" t '|a|_1, '.ftitle(a11, b11), \
     exp (f11(log(x))) ls 1 lc rgb "coral" notitle, \
     exp (f111(log(x))) ls 1 lc rgb "pink" t ftitle(a111, b111), \
     '../tourniquet/log' u 1:3 w p pt 7 ps 1.5 lc rgb "navy" t '|a|_2, '.ftitle(a2, b2), \
     exp (f2(log(x))) ls 1 lc rgb "red" notitle, \
     'log' u (2**$1):3 w p pt 7 ps 1.5 lc rgb "dark-green" t '|a|_2, '.ftitle(a22, b22), \
     exp (f222(log(x))) ls 1 lc rgb "pink" t ftitle(a222, b222), \
     exp (f22(log(x))) ls 1 lc rgb "coral" notitle, \
     '../tourniquet/log' u 1:4 w p pt 5 ps 1.5 lc rgb "skyblue" t '|a|_{max}, '.ftitle(a3, b3), \
     exp (f3(log(x))) ls 1 lc rgb "red" notitle, \
     'log' u (2**$1):4 w p pt 5 ps 1.5 lc rgb "forest-green" t '|a|_{max}, '.ftitle(a33, b33), \
     exp (f33(log(x))) ls 1 lc rgb "coral" notitle

Spatial convergence for a (script)

reset
set xlabel 'N_{max}'
set ylabel 'L_1(q),L_2(q),L_{max}(q)'
set format y '%.1e'
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) '../tourniquet/log' u (log($1)):(log($5)) via a1, b1
fit f2(x) '../tourniquet/log' u (log($1)):(log($6)) via a2, b2
fit f3(x) '../tourniquet/log' u (log($1)):(log($7)) via a3, b3

f11(x) = a11 + b11*x
f22(x) = a22 + b22*x
f33(x) = a33 + b33*x
fit f11(x) 'log' u (log(2**$1)):(log($5)) via a11, b11
fit f22(x) 'log' u (log(2**$1)):(log($6)) via a22, b22
fit f33(x) 'log' u (log(2**$1)):(log($7)) via a33, b33

f111(x) = a111 + b111*x
f222(x) = a222 + b222*x
fit [x=2^10:] f111(x) 'log' u (log(2**$1)):(log($5)) via a111, b111
fit [x=2^10:] f222(x) 'log' u (log(2**$1)):(log($6)) via a222, b222

plot '../tourniquet/log' u 1:5 w p pt 6 ps 1.5 lc rgb "blue" t '|q|_1, '.ftitle(a1, b1), \
     exp (f1(log(x))) ls 1 lc rgb "red" notitle, \
     'log' u (2**$1):5 w p pt 6 ps 1.5 lc rgb "sea-green" t '|q|_1, '.ftitle(a11, b11), \
     exp (f11(log(x))) ls 1 lc rgb "coral" notitle, \
     exp (f111(log(x))) ls 1 lc rgb "pink" t ftitle(a111, b111), \
     '../tourniquet/log' u 1:6 w p pt 7 ps 1.5 lc rgb "navy" t '|q|_2, '.ftitle(a2, b2), \
     exp (f2(log(x))) ls 1 lc rgb "red" notitle, \
     'log' u (2**$1):6 w p pt 7 ps 1.5 lc rgb "dark-green" t '|q|_2, '.ftitle(a22, b22), \
     exp (f22(log(x))) ls 1 lc rgb "coral" notitle, \
     exp (f222(log(x))) ls 1 lc rgb "pink" t ftitle(a222, b222), \
     '../tourniquet/log' u 1:7 w p pt 5 ps 1.5 lc rgb "skyblue" t '|q|_{max}, '.ftitle(a3, b3), \
     exp (f3(log(x))) ls 1 lc rgb "red" notitle, \
     'log' u (2**$1):7 w p pt 5 ps 1.5 lc rgb "forest-green" t '|q|_{max}, '.ftitle(a33, b33), \
     exp (f33(log(x))) ls 1 lc rgb "coral" notitle

Spatial convergence for q (script)