A Flow with Shocks :

Compressible Euler equations should be discretized with Finite Volumes or FEM with flux up-winding scheme but these are not implemented in FreeFem++. Nevertheless acceptable results can be obtained with the method of characteristics provided that the mean values \(\bar f=\frac12(f^++f^-)\) are used at shocks in the scheme, and finally mesh adaptation .
\( \partial_t\rho+\bar u\nabla\rho + \bar\rho\nabla\cdot u=0 \)
\( \bar\rho( \partial_t u+\frac{\overline{\rho u}}{\bar\rho}\nabla u +\nabla p=0 \)
\( \partial_t p + \bar u\nabla p +(\gamma-1)\bar p\nabla\cdot u =0 \)

[euler] One possibility is to couple \(u,p\) and then update \(\rho\), i.e.

\( \frac 1{(\gamma-1)\delta t\bar p^m} (p^{m+1}-p^m \circ X^m) + \nabla\cdot u^{m+1} =0 \)
\( \frac{\bar\rho^m}{\delta t}(u^{m+1}-u^m \circ {\tilde X}^m ) +\nabla p^{m+1}=0 \)
\( \rho^{m+1} = \rho^m \circ X^m + \frac{\bar\rho^m}{(\gamma-1)\bar p^m}(p^{m+1}-p^m \circ X^m) \)

[eulalgo] A numerical result is given on Figure [figvfive] and the FreeFem++ script is


 [figvfive] Pressure for a Euler flow around a disk at Mach 2 computed by ([eulalgo]) 

 int anew=1;
 real x0=0.5,y0=0, rr=0.2;
 border ccc(t=0,2){x=2-t;y=1;};
 border ddd(t=0,1){x=0;y=1-t;};
 border aaa1(t=0,x0-rr){x=t;y=0;};
 border cercle(t=pi,0){ x=x0+rr*cos(t);y=y0+rr*sin(t);}
 border aaa2(t=x0+rr,2){x=t;y=0;};
 border bbb(t=0,1){x=2;y=t;};

 int m=5; mesh Th;
 if(anew) Th = buildmesh (ccc(5*m) +ddd(3*m) + aaa1(2*m) + cercle(5*m)
               + aaa2(5*m) + bbb(2*m) );
       else Th = readmesh("Th_circle.mesh"); plot(Th,wait=0);

 real dt=0.01, u0=2, err0=0.00625, pena=2;
 fespace Wh(Th,P1);
 fespace Vh(Th,P1);
 Wh u,v,u1,v1,uh,vh;
 Vh r,rh,r1;
 macro dn(u) (N.x*dx(u)+N.y*dy(u) ) //  def the normal derivative

 if(anew){ u1= u0; v1= 0; r1 = 1;}
 else {
     ifstream g("u.txt");g>>u1[];
     ifstream gg("v.txt");gg>>v1[];
     ifstream ggg("r.txt");ggg>>r1[];
     plot(u1,ps="eta.eps", value=1,wait=1);
     err0=err0/10; dt = dt/10;

 problem  eul(u,v,r,uh,vh,rh)
    = int2d(Th)(  (u*uh+v*vh+r*rh)/dt
                   + ((dx(r)*uh+ dy(r)*vh) - (dx(rh)*u + dy(rh)*v))
  + int2d(Th)(-(rh*convect([u1,v1],-dt,r1) + uh*convect([u1,v1],-dt,u1)
                 + vh*convect([u1,v1],-dt,v1))/dt)
   +int1d(Th,6)(rh*u)   // +int1d(Th,1)(rh*v)
  + on(2,r=0) + on(2,u=u0) + on(2,v=0);

 int j=80;
 for(int k=0;k<3;k++)
     if(k==20){ err0=err0/10; dt = dt/10; j=5;}
     for(int i=0;i<j;i++){
        eul; u1=u; v1=v; r1=abs(r);
         cout<<"k="<<k<<"  E="<<int2d(Th)(u^2+v^2+r)<<endl;
 Th = adaptmesh (Th,r, nbvx=40000,err=err0,
       abserror=1,nbjacoby=2, omega=1.8,ratio=1.8, nbsmooth=3,
       splitpbedge=1, maxsubdiv=5,rescaling=1) ;

 ofstream f("u.txt");f<<u[];
 ofstream ff("v.txt");ff<<v[];
 ofstream fff("r.txt");fff<<r[];
 r1 = sqrt(u*u+v*v);
 plot(r1,ps="mach.eps", value=1);