//  two_fluid.edp
//
//  Discussion:
//
//    Two−fluid incompressible flow
//
//    div u = 0 in Omega
//    D_t rho = 0 in Omega
//    D_t(rho u) + div (mu grad u) + grad p = rho g.
//
//    MU is the dynamic viscosity,
//    MU = RHO * NU, where NU is the kinematic viscosity.
//
//    D_t is the material derivative:
//
//      D_t rho = d/dt rho + d/dx rho * v
//
//    where v is the convection velocity.
//
//    The fluids have different densities, 1 and 10, so the local
//    density, between 1 and 10, reflects the relative proportion of
//    the two fluids.
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/two_fluid/two_fluid.edp
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    25 June 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
//
cout << "\n";
cout << "two_fluid\n";
cout << "  FreeFem++ version\n";
cout << "  Flow of two incompressible fluids of different densitites.\n";

real gravity = 9.81;
real rhom = 1.0;
real rhop = 10.0;
real num = 0.05;
real nup = 0.01;
real dt = 0.2;
//
//  Define a fine mesh.
//
real Lx = 2;
real Ly = 1;
mesh Th = square ( 120, 60, [ x * Lx, y * Ly ] );
plot ( Th, cmm = "two_fluid Mesh", wait = true, ps = "two_fluid_mesh.ps" );
//
//  Define the finite element space.
//
fespace Uh ( Th, P1b );
fespace Vh ( Th, P1 );
fespace Wh ( Th, P2 );
//
//  Define test and trial variables.
//
Uh u1;
Uh u2;
Uh u1old;
Uh u2old;
Uh u1h;
Uh u2h;
Vh p;
Vh ph;
Wh phi;
Wh phih;
Wh phiold;
Wh rho;
Wh rhoold;
Wh mu;
//
// Initialization
//
phi = 0.75 * Ly / Lx * x - y;
phiold = phi;
u1 = 0.0;
u1old = u1;
u2 = 0.0;
u2old = u2;
//
//  Define the weak form of the convection equation.
//
problem ConvectLevel ( phi, phih ) =
    int2d(Th) ( phi * phih / dt )
  - int2d(Th) ( convect ( [u1old,u2old], -dt, phiold ) * phih / dt )
  + int2d(Th) ( 0.001 * dx(phi)*dx(phih) + 0.001 * dy(phi)*dy(phih) );
//
//  Define the weak form of the Navier Stokes equations.
//
problem NavierStokes ( [u1, u2, p], [u1h, u2h, ph] ) =
    int2d(Th) (rho*u1*u1h / dt)
  - int2d(Th) (rho*convect([u1old,u2old], -dt, u1old)*u1h/dt)
  + int2d(Th) (mu*dx(u1)*dx(u1h)+mu*dy(u1)*dy(u1h))
  - int2d(Th) (p*dx(u1h))
  + int2d(Th) (rho*u2*u2h/dt)
  - int2d(Th) (rho*convect([u1old,u2old], -dt, u2old)*u2h/dt)
  + int2d(Th) (mu*dx(u2)*dx(u2h)+mu*dy(u2)*dy(u2h))
  - int2d(Th) (p*dy(u2h))
  + int2d(Th) (rho*gravity*u2h)
  + int2d(Th) ((dx(u1)+dy(u2))*ph+0.000001*p*ph)
  + on ( 1, 2, 3, 4, u1=0.0, u2=0.0 );
//
for ( int it = 0; it < 30; it++ )
{
  for ( int innerloop = 0; innerloop < 1; innerloop++ )
  {
//
//  Convect PHI using the current velocity field.
//
    phiold = phi;
    ConvectLevel;
//
//  Update RHO and MU based on the value of PHI.
//
    rhoold = rho;
    rho = rhom * (phi<=0) + rhop * (phi>0);

    mu = rhom * num * (phi<=0) + rhop * nup * (phi>0);
//
//  Solve the Navier Stokes equations for the new velocity field.
//
    NavierStokes;
    u1old = u1;
    u2old = u2;
  }
//
//  Plot the density field.
//
  plot ( rho, nbiso = 40, fill = true, value = true, 
    cmm = "two_fluid iteration " + it,
    ps = "two_fluid_rho_"+it+".ps" );

  cout << "  IT = " << it 
       << ", pmin = " << p[].min 
       << ", pmax = " << p[].max << endl;
}
//
//  Terminate.
//
cout << "\n";
cout << "two_fluid:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );