/**
# 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:
$$
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:
$$
\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:
$$
\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);
}

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

/**
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 different grid sizes. */

  for (N = 128; N <= 1024; N *= 2) {
    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 (N == 1024) {

    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\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[]
	       );
  }
}

/**
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",
	   N,
	   na.avg, na.rms, na.max,
	   nq.avg, nq.rms, nq.max);
}

/**
## End of simulation */

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

/**
# Results for first order

#### Cross-sectional area and flow rate

We first plot the spatial evolution of the cross-sectional area $a$ at
$t={0, 0.01, 0.02, 0.03, 0.04}$ for $N=1024$.

~~~gnuplot $a/a_0$ for $N=1024$. 
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'
~~~

~~~gnuplot $q$ for $N=1024$. 
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'
~~~

#### 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$.

~~~gnuplot Spatial convergence for $a$
reset
set xlabel 'N'
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) 'log' u (log($1)):(log($2)) via a1, b1
fit f2(x) 'log' u (log($1)):(log($3)) via a2, b2
fit f3(x) 'log' u (log($1)):(log($4)) via a3, b3

plot '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 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 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
~~~

~~~gnuplot Spatial convergence for $q$
reset
set xlabel 'N'
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) 'log' u (log($1)):(log($5)) via a1, b1
fit f2(x) 'log' u (log($1)):(log($6)) via a2, b2
fit f3(x) 'log' u (log($1)):(log($7)) via a3, b3

plot '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 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 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
~~~
*/