//  shock.edp
//
//  Discussion:
//
//    Supersonic perfect gas flow around an ellipse.
//
//    Be careful : nonconservative formulation of the compressible 
//    Euler equations.
//
//    A conservative Finite Volume method should be used.
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/shock/shock.edp
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    20 May 2015
//
//  Author:
//
//    Florian De Vuyst
//
//  Reference:
//
//    Florian De Vuyst,
//    Numerical modeling of transport problems using freefem++ software -
//    with examples in biology, CFD, traffic flow and energy transfer,
//    HAL id: cel-00842234
//    https://cel.archives-ouvertes.fr/cel-00842234
//
//  Parameters:
//
//    real ALPHA, the angle of attack, in radians.
//
//    real CINF, the speed of sound.
//
//    real DT, time step.
//
//    real RA, major radius of the ellipse.
//
//    real RADIUS, the radius of the infinite flow boundary.
//
//    real RB, minor radius of the ellipse.
//
//    real UINF, the infinite flow Mach number is 2.
//
//    real X0, X translation of the body.
//
cout << "\n";
cout << "shock\n";
cout << "  FreeFem++ version\n";
cout << "  Solve the Euler equations for perfect gas flow\n";
cout << "  around an elliptical obstacle.\n";

real gamma = 1.4;
real pinf = 100.0;
real rhoinf = 0.3;
real cinf = sqrt ( gamma * pinf / rhoinf );
real uinf = 1.5 * cinf;
real radius = 10.0;
real alpha = 0.3;
real x0 = -4.0;
real ra = 2.0;
real rb = 0.3;
real dt = 0.01;
//
//  External "infinite" boundary is a circle
//
border binf (t=0,2*pi) {x=radius*cos(t); y=radius*sin(t);}
//
//  Obstacle is an ellipse
//
border wall (t=0,2*pi) { x = x0+cos(alpha)*ra*cos(t)+sin(alpha)*rb*sin(t);
                         y = -sin(alpha)*ra*cos(t)+cos(alpha)*rb*sin(t);}
//
//  Define the mesh.
//
mesh Th = buildmesh ( binf ( 100 ) + wall ( -80 ) );
plot ( Th, wait = true, ps = "shock_mesh.ps" );
//
//  Define the finite element space.
//
fespace Vh ( Th, P1 );
Vh ap;
Vh arho;
Vh p;
Vh rho;
Vh T;
Vh u1;
Vh u2;
Vh rhoold, pold, Told, arhoold, apold, u1old, u2old;
//
//  Define test functions.
//
Vh v;
Vh v1;
Vh v2;
Vh w;
//
//  Field initialization
//
rho = rhoinf; 
p = pinf;
u1 = uinf;
u2 = 0;
arho = log ( rho );
ap = log ( p );
T = p / rho; // temperature−like

plot ( T, nbiso = 50, fill = true, value = true, wait = false );
//
//  Define the weak form of the Euler equations.
//
problem euler ( [arho, u1, u2, ap], [v, v1, v2, w] ) =
  int2d ( Th ) ( arho*v/dt )
- int2d ( Th ) ( convect([u1old,u2old], -dt, arhoold)*v/dt )
+ int1d ( Th, binf) ( (uinf*N.x)*v )
- int2d ( Th ) ( u1*dx(v)+u2*dy(v) )
+ int2d ( Th ) ( u1*v1/Told/dt )
- int2d ( Th ) ( convect([u1old,u2old], -dt, u1old)*v1/Told/dt )
+ int2d ( Th ) ( u2*v2/Told/dt )
- int2d ( Th ) ( convect([u1old,u2old], -dt, u2old)*v2/Told/dt )
+ int2d ( Th ) ( dx(ap)*v1+dy(ap)*v2 )
+ int2d ( Th ) ( ap*w/dt )
- int2d ( Th ) ( convect([u1old,u2old], -dt, apold)*w/dt )
+ int1d ( Th, binf) ( gamma*(uinf*N.x)*w )
- int2d ( Th ) ( gamma*u1*dx(w)+gamma*u2*dy(w) );

for ( int it = 0; it < 80; it++ )
{
  u1old = u1;
  u2old = u2;
  arhoold = arho;
  apold = ap;
  Told = T;
  euler;
  if ( 70 < it )
  {
    Th = adaptmesh ( Th, ap );
    u1 = u1;
    u2 = u2;
    arho = arho;
    ap = ap;
  }
  rho = exp ( arho );
  p = exp ( ap );
  T = p / rho;
  plot ( rho, nbiso = 50, fill = false, value = true );
  if ( it < 20 ) 
  {
    plot ( rho, nbiso = 50, fill = true, value = true );
  }
}
//
//  Plot final results.
//
plot ( T, nbiso = 50, fill = true, value = true, ps = "shock_t.ps" );
plot ( T, nbiso = 50, fill = false, value = true, ps = "shock_tiso.ps" );
plot ( p, nbiso = 50, fill = true, value = true, ps = "shock_p.ps" );
plot ( rho, nbiso = 50, fill = true, value = true, ps = "shock_rho.ps" );
plot ( [ u1, u2 ], ps = "shock_velocity.ps" );
plot ( Th, ps = "shock_mesh2.ps" );
//
//  Terminate.
//
cout << "\n";
cout << "shock:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );