sandbox/lopez/unsymshock.c
Unsymmetrical explosion
We solve the Euler equations for a compressible gas under a nonsymmetric initial pressure (density) distribution. This problem has been recently aborded by Eggers et al. (2017). These authors showed that a selfsimilar structure exists in the vicinity of the emerging shock wave. The shock spreads in the transversal direction as |t_o −t|^{1/2} and along the direction of propagation as |t_o − t|^{3/2}. t_o is the singularity time.
#include "compressible.h"
#define LEVEL 8
int main() {
We set Neumann conditions for normal velocity at all the boundaries.
foreach_dimension() {
w.n[right] = neumann(0);
w.n[left] = neumann(0);
}
The domain spans [-2:2]\times[-2:2].
origin (-2, -2);
size (4.);
init_grid (1 << (LEVEL));
run();
}
The initial density distribution is \displaystyle \rho(x,y,0) = 0.2 + e^{-4(x^4 + y^2)}
The momentum, \mathbf{w} = \rho \mathbf{u}, is initially null while the total energy is E = \rho \mathbf{u}^2/2 + p/(\gamma-1). We assume that, initially, the isentropic relationship, p/\rho^\gamma = cte, holds.
foreach() {
rho[] = 0.2 + exp(-4*(sq(y) + sq(x)*sq(x)));
foreach_dimension()
w.x[] = 0.;
E[] = pow(rho[], gammao)/(gammao - 1.);
}
}
We wish to plot density distributions at different times.
event plot (t = {0., 0.4, 0.511, 0.55}) {
char name[80];
sprintf (name, "density_%g", t);
FILE * fp = fopen (name, "w");
output_matrix(rho, fp, 1000);
fclose (fp);
}
We will plot the time evolution of maximum entropy and the inverse of the maximum gradient of the density.
event graphs (i++) {
scalar rhograd[]; //gradient modulus
scalar S[]; //Entropy
foreach() {
double a = 0., w2 = 0.;
foreach_dimension() {
a += sq(rho[1] -rho[-1]);
w2 += sq(w.x[]);
}
rhograd[] = sqrt(a)/(2.*Delta);
double pres = (E[]-0.5*w2/rho[])*(gammao-1.); //pressure
S[] = log(pres/(pow(rho[],(gammao))));
}
stats sr = statsf (rhograd);
stats ss = statsf (S);
fprintf (stderr, "%g %g %g\n", t, 1./sr.max, ss.max);
}
We adapt the mesh by controlling the error on the density field.
Results
We reproduce here figure 2 and 5 (sublots a and b) of the works of Eggers et al. (2017). Note that figures of the paper were computed with a largely finer mesh.
set pm3d
#set hidden
set cntrparam levels 7
set view 55,40
unset key
unset map
unset colorbox
set palette color
unset tics
set xyplane relative -1.2
set contours
set multiplot layout 2,2
splot [0:2.][-1:1][0:1] 'density_0' binary u 1:2:3 with pm3d
splot [0:2.][-1:1][0:1] 'density_0.4' binary u 1:2:3 with pm3d
splot [0:2.][-1:1][0:1] 'density_0.511' binary u 1:2:3 with pm3d
splot [0:2.][-1:1][0:1] 'density_0.55' binary u 1:2:3 with pm3d
unset multiplot
set key left
set multiplot layout 2,1
set ylabel '|grad^{-1} ({/Symbol r})|_{max}'
set xtics 0.02
set ytics 0.5
plot [0.4:0.56][0:0.6] 'log' u 1:2 w l t 'Basilisk', 3.85*(0.511-x) t '|t - t_o| law'
set ylabel 'S_{max}'
set xlabel 't'
set ytics 0.005
plot [0.4:0.56][0:0.009] 'log' u 1:3 w l t 'Basilisk', 0.7*(x-0.511)**1.5 t '|t - t_o|^{3/2} law'
unset multiplot