//  thermal_convection.edp
//
//  Discussion:
//
//    Forced + Natural heat convection in a pipe
//    Navier−Stokes equations + convection−diffusion on temperature
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/thermal_convection/thermal_convection.edp
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    26 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 << "thermal_convection:\n";
cout << "  FreeFem++ version\n";
cout << "  Thermal convection and diffusion by Navier-Stokes flow.\n";
//
//  Define the boundary lines.
//
real lx = 0.25;
real Lx = 3.0;
real Ly = 1.0;

border c1 (t=0.0, lx)   {x = t;   y = 0.0;}
border c2 (t=lx, Lx-lx) {x = t;   y = 0.0;}
border c3 (t=Lx-lx,Lx)  {x = t;   y = 0.0;}
border c4 (t=0.0,Ly)    {x = Lx;  y = t;}
border c5 (t=Lx,0.0)    {x = t;   y = Ly;}
border c6 (t=Ly,0.0)    {x = 0.0; y = t;}
//
//  Define the mesh.
//
mesh Th = buildmesh ( c1(10) + c2(200) + c3(10) + c4(30) + c5(100) + c6(30) );
//
//  Plot the mesh.
//
plot ( Th, wait = false, ps = "thermal_convection_mesh.ps" );
//
//  Define the finite element spaces:
//  UH for velocity fields,
//  XH for pressure P, temperature THETA, and specific volume TAU.
//
fespace Uh(Th, P1b );
Uh u, v, uold, vold, uh, vh;
fespace Xh(Th, P1 );
Xh p, ph;
Xh theta, thetaold, thetah;
Xh tau;
//
//  Initialize the fields.
//
real uin0 = 0.2;
func uin = uin0 * 4.0 * (y/Ly) * (1.0-y/Ly);
u = uin;
uold = u;
v = 0.0;
vold = v;

real thetain = 20.0;
theta = thetain;
thetaold = theta;
//
//  Define the weak form of the heat equation.
//
real dt = 0.05;
real heatflux = 100.0;
real kappa = 0.001;

problem HeatStep ( theta, thetah ) =
    int2d (Th) ( theta * thetah / dt )
  - int2d (Th) ( convect ([u,v], -dt, thetaold)*thetah/dt )
  + int2d (Th) ( kappa * dx(theta) * dx(thetah) + kappa * dy(theta) * dy(thetah) )
  - int1d (Th, c2) ( kappa * heatflux * thetah )
  + on ( c6, theta = thetain );
//
//  Define the weak form of the Navier Stokes equations.
//
real gravity = 9.81;
real nu = 0.001;

problem NavierStokesStep ( [u, v, p], [uh, vh, ph] ) =
    int2d (Th) ( u * uh / dt )
  - int2d (Th) ( convect ( [uold, vold], -dt, uold ) * uh/dt )
  + int2d (Th) ( nu * dx(u) * dx(uh) + nu * dy(u) * dy(uh) )
  + int2d (Th) ( tau * dx(p) * uh )
  + int2d (Th) ( v * vh / dt )
  - int2d (Th) ( convect ([uold, vold], -dt, vold) * vh / dt )
  + int2d (Th) ( nu * dx(v) * dx(vh) + nu * dy(v) * dy(vh) )
  + int2d (Th) ( tau * dy(p) * vh )
  + int2d (Th) ( gravity * vh )
  + int2d (Th) ( dx(u) * ph + dy(v) * ph )
  + on ( c6, u=uin, v=0 )
  + on ( c1, c2, c3, c5, u=0, v=0 );
//
//  Time steps
//
real cst = 1.0;
real dtsnap = 0.5;
real t = 0.0;
real ttosnap = 2.0;

for ( int it = 1; it <= 100; it++ )
{
  t = it * dt;
  cout << it << "  T = " << t << "\n";

  HeatStep;

  thetaold = theta;
  tau = cst * theta;

  NavierStokesStep;

  uold = u;
  vold = v;

  plot ( [ u, v ], theta, nbiso = 80, value = false );

  if ( ttosnap <= t )
  {
    cout << "  Plotting T = " << t << "\n";
    ttosnap = ttosnap + dtsnap;
    plot ( [u,v], theta, nbiso = 80, value = false, 
      ps = "thermal_convection_"+it+".ps",
      cmm = "Step "+it );
  }
}
//
//  Terminate.
//
cout << "\n";
cout << "thermal_convection:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );