/**
# Two- and three-dimensional explosions

We solve the Euler equations for a compressible gas. We also need to
compute volume fractions for the initial condition. */

#include "compressible.h"
#include "fractions.h"

#if dimension == 2
# define LEVEL 7
#else // 3D
# define LEVEL 6
#endif

int main() {

  /**
  We make boundary conditions free outflow. */

  foreach_dimension() {
    w.n[right] = neumann(0);
    w.n[left]  = neumann(0);
  }
  
  /**
  The domain spans $[-1:1]\times[-1:1]\times[-1:1]$. */

  origin (-1, -1, -1);
  size (2.);
  init_grid (1 << LEVEL);
  run(); 
}

/**
Initial conditions come from Toro's book (Riemann Solvers and
Numerical Methods for Fluid Dynamics, 3rd Edition, Springer Ed.)
Chapter 17 section 17.1.1 are given in terms of density ($\rho$),
pression ($p$), velocity ($u$) both at the left and right side of the
discontinuity placed at $R=0.4$. */

event init (t = 0)
{
  double R = 0.4 ;
  double rhoL = 1., rhoR = 0.125 ;
  double pL = 1.0,  pR = 0.1 ;
  
  /**
  To define an accurate (i.e. symmetrical) initial sphere of rayon
  $R$, we compute the volume fraction corresponding to the spherical
  interface. */

  scalar f[];
  fraction (f, sq(x) + sq(y) + sq(z) - sq(R));
  
  /**
  Left and right initial states for $\rho$, $\mathbf{w}$ and energy
  $E = \rho \mathbf{u}^2/2 + p/(\gamma-1)$. */
  
  foreach() {
    rho[] = rhoR*f[] + rhoL*(1. - f[]);
    foreach_dimension()
      w.x[] = 0.;
    E[] = (pR*f[] + pL*(1. - f[]))/(gammao - 1.);
  }

  theta = 1.3; // tune limiting from the default minmod
}

event print (t = 0.25)
{

  /**
  At $t=0.25$ we output the values of $\rho$ and the normal velocity
  $\mathbf{u}_n$ as functions of the radial coordinate. */

  foreach() {
    double r = sqrt(sq(x) + sq(y) + sq(z));
    double wn = (w.x[]*x + w.y[]*y + w.z[]*z)/r;
    printf ("%g %g %g\n", r, rho[], wn/rho[]);
  }

  /**
  For reference we also output a cross-section at $y=0$. */

  for (double x = 0; x <= 1; x += 1e-2)
    fprintf (stderr, "%g %.4f %.4f\n", x,
	     interpolate (rho, x, 0.),
	     interpolate (w.x, x, 0.));
}

/**
On trees, we adapt the mesh by controlling the error on the density
field. */

#if TREE
event adapt (i++) {
  adapt_wavelet ({rho}, (double[]){5e-3}, LEVEL + 1);
}
#endif

/**
## Results

Results are presented in terms of $\rho$ and normal velocity $u_n$ for
Cartesian (7 levels in 2D and 6 levels in 3D) and adaptive (8 levels
in 2D and 7 levels in 3D) computations. The numerical results compare
very well with Toro's numerical experiments.

~~~gnuplot Radial density profile
set xrange [0:1]
set xlabel 'r'

set term PNG enhanced font ",10"
set output 'rho.png'
set ylabel 'rho'
plot './cout' u 1:2 w p pt 7 ps 0.2 t '2D Cartesian', \
     './out' u 1:2 w p pt 7 ps 0.2 t '2D Adaptive', \
     '../explosion3D/out' u 1:2 w p pt 7 ps 0.2 t '3D Cartesian', \
     '../explosion.3D/out' u 1:2 w p pt 7 ps 0.2 t '3D Adaptive'
~~~

~~~gnuplot Normal velocity
set output 'velocity.png'
set ylabel 'Normal velocity'
plot './cout' u 1:3 w p pt 7 ps 0.2 t '2D Cartesian',		  \
     './out' u 1:3 w p pt 7 ps 0.2 t '2D Adaptive',		  \
     '../explosion3D/out' u 1:3 w p pt 7 ps 0.2 t '3D Cartesian', \
     '../explosion.3D/out' u 1:3 w p pt 7 ps 0.2 t '3D Adaptive'
~~~

*/