/**
# Advection of two fluids at different pressures

This problem, initially proposed by [Johnsen and Colonius,
2006](#johnsen2006) tests the solver for the advection of an interface
between two different ideal gases at uniform velocity and pressure
quantifying the amplitude of spurious pressure and velocity
oscillations induced by the method when advecting an interface with
different material properties. See also Figure 4 in [Fuster and
Popinet, 2018](#fuster2018).

$$(\rho,u,p,\gamma)^T_L = (1, 0.5, 1/1.4, 1.2)^T$$
$$(\rho,u,p,\gamma)^T_R = (10, 0.5, 1/1.4, 1.4)^T$$

By construction the method should keep the pressure and velocity
uniform. During the advection step the energy is advected avoiding any
diffusion and therefore both the provisional pressure and velocity are
uniform and do not need to be corrected during the projection step. */

#include "grid/multigrid1D.h"
#include "compressible/two-phase.h"
#include "compressible/Mie-Gruneisen.h"

/**
We set the problem parameters, size domain and boundary conditions. */

double rhoL = 1., rhoR = 10.;
double pL, pR;
double tend = 8;

int main()
{
  L0 = 2. [1];
  X0 = - L0/2.;
  
  periodic (right);
  
  pL = pR = 1./1.4;
  gamma1 = 1.2;
  gamma2 = 1.4;
 
  N = 128;
  run();
}

/** 
The initial conditions are: */

event init (i = 0)
{  
  double u0 = 0.5;
  foreach() {
    f[] = (x < 0.);
    p[] = f[]*pL + (1. - f[])*pR;
    frho1[] = f[]*rhoL;
    frho2[] = (1. - f[])*rhoR;
    fE1[] = f[]*pL/(gamma1 - 1.) +  0.5*frho1[]*sq(u0);
    fE2[] = (1. - f[])*pR/(gamma2 - 1.) +  0.5*frho2[]*sq(u0);
    q.x[] = (frho1[] + frho2[])*u0;
  }
}

/** 
We output the field variables at the end of the simulation. */

event outputdata (t = tend)
{  
  scalar perr[], uerr[];
  
  foreach () {
    perr[] = fabs(p[] - 1./1.4);
    uerr[] = fabs(q.x[]/rho[] - 0.5);

    double Ek = 0.;
    foreach_dimension()
      Ek += sq(q.x[]);

    fprintf (stderr, "%g %g %g %g %g %g %g \n", x, t, p[], rho[], q.x[]/rho[], f[],
	     (fE1[] + fE2[] - 0.5*Ek/rho[])/(f[]/(gamma1 - 1.) + (1. - f[])/(gamma2 - 1.)));
  }

  stats sp = statsf(perr), su = statsf(uerr);
  fprintf (stdout, "%g %g\n", sp.sum/sp.volume, su.sum/su.volume);
  fflush (stdout);
  assert (sp.sum/sp.volume < 2e-9 && su.sum/su.volume < 2e-9);
}

/**
The results below are those displayed in Figure 4 of [Fuster and
Popinet, 2018](#fuster2018).

~~~gnuplot Pressure profile
set xlabel 'x'
set ylabel 'p'
unset key
p "log" u 1:3 w lp
~~~ 
 
~~~gnuplot Density profile
set xlabel 'x'
set ylabel '{/Symbol r}'
p "log" u 1:4 w lp
~~~ 

~~~gnuplot Velocity profile
set xlabel 'x'
set ylabel 'u'
p "log" u 1:5 w lp
~~~ 

~~~gnuplot VOF profile
set xlabel 'x'
set ylabel 'VOF function'
p "log" u 1:6 w lp
~~~

### References

~~~bib
@article{johnsen2006,
title = {Implementation of {WENO} schemes in compressible multicomponent flow problems},
journal = {Journal of Computational Physics},
volume = {219},
number = {2},
pages = {715-732},
year = {2006},
issn = {0021-9991},
doi = {https://doi.org/10.1016/j.jcp.2006.04.018},
url = {https://www.sciencedirect.com/science/article/pii/S0021999106002014},
author = {Eric Johnsen and Tim Colonius}
}

@hal{fuster2018, hal-01845218}
~~~
*/