/**
# Inviscid wave propagation in a straight elastic artery

We solve the 1D blood flow equations in a straight artery initially
deformed in its center. At $t>0$, the vessel relaxes towards its
steady-state at rest. Consequently, two waves are created that
propagate towards both extremeties of the artery. The solution of the
flow of blood generated by this elastic relaxation is obtained
numerically and compared to the anaytic solution obtained using linear
wave theory. */

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

/**
We define the artery's geometrical and mechanical properties. */

#define R0 (1.)
#define DR (1.e-3)
#define XS (2./5.*(L0))
#define XE (3./5.*(L0))
#define shape(x) (((XS) <= (x) && (x) <= (XE)) ? (DR)/2.*(1. + cos(pi + 2.*pi*((x) - (XS))/((XE) - (XS)))) : 0.)

#define K0 1.e4
#define C0 (sqrt (0.5*(K0)*sqrt(pi)*(R0)))

/**
We define the linear analytic solution for the cross-sectional area
$a$ and the flow rate $q$, which depends on the shape of initial
vessel deformation. */

#define analytic_a(t,x) (pi*sq (R0)*sq (1. + 0.5*((shape ((x) - (C0)*(t))) + (shape ((x) + (C0)*(t))))))
#define analytic_q(t,x) (-(C0)*(-(shape ((x) - (C0)*(t))) + (shape ((x) + (C0)*(t))))*(analytic_a((t),(x))))

int main() {

  /**
  The domain is 10.*/

  L0 = 10.;
  size (L0);
  origin (0.);
  
  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[]*(1. + (shape (x))));
    q[] = 0.;
  }
}

/**
## Post-processing

We output the computed fields. */
 
event output (t = {0., 0.01, 0.02, 0.03, 0.04}) {

  if (N == 256) {
    
    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))/(pi*sq ((R0))), a[]/(pi*sq ((R0))),
	       (analytic_q (t,x)), q[]
	       );
  }
}

/**
Next, we compute the spatial error for the flow rate. */

event error (t = 0.03) {

  scalar err_q[];
  foreach()
    err_q[] = fabs (q[] - (analytic_q (t, x)));
  boundary ((scalar *) {err_q});

  norm nq = normf (err_q);
  
  fprintf (ferr, "%d %g %g %g\n",
	   N,
	   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=256$.

~~~gnuplot $a/a_0$ for $N=256$. 
reset
set xlabel 'x'
set ylabel 'a/a_0'
set yrange [1:1 + 2.5e-3]
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=256$. 
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 flow rate $q$ with
the number of cells $N$.

~~~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($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 '|q|_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 '|q|_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 '|q|_{max}, '.ftitle(a3, b3), \
     exp (f3(log(x))) ls 1 lc rgb "red" notitle
~~~
*/