//  Discussion:
//
//    Create a movie of the heat transfer in a von Karman vortex street.
//
//    Export the data, using ffmatlib, for graphic processing by MATLAB or Octave.
//
//  Licensing:
//
//    Copyright (C) 2018 Chloros2 <chloros2@gmx.de>
//
//    This program is free software: you can redistribute it and/or modify it
//    under the terms of the GNU General Public License as published by
//    the Free Software Foundation, either version 3 of the License, or
//    (at your option) any later version.
//
//    This program is distributed in the hope that it will be useful, but
//    WITHOUT ANY WARRANTY; without even the implied warranty of
//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//    GNU General Public License for more details.
//
//    You should have received a copy of the GNU General Public License
//    along with this program.  If not, see
//    <https://www.gnu.org/licenses/>.
//
//  Modified:
//
//    2018-05-19
//
// Author: 
//
//   Chloros2 <chloros2@gmx.de>
//
include "ffmatlib.idp"

cout << endl;
cout << "karman_vortex:" << endl;
cout << "  FreeFem++ version." << endl;
cout << "  Compute fluid flow past an obstacle." << endl;
cout << "  Observe heat transfer in the resulting von Karman vortex street." << endl;
cout << "  Using the ffmatlib interface, export the data for graphic" << endl;
cout << "  processing by MATLAB/Octave." << endl;

real scale = 0.5;
real width = 4.0 * scale;
real height = 0.8 * scale;
real R = 0.04 * scale;

int C1 = 1;
int C2 = 2;
int C3 = 3;
int C4 = 4;
int C5 = 5;
int C6 = 6;

border floor   ( t = 0.0, width )
  { x = t; y = 0.0; label = C6; };
border ceiling ( t = width, 0.0 )
  { x = t; y = height; label = C1; };
border right   ( t = 0, height )
  { x = width; y = t; label = C3; };
border left    ( t = height, 0.0 )
  { x = 0.0; y = t; label = C2; };
border cir     ( t = 2.0 * pi, 0.0 )
  { x = 0.1 * width + R * cos ( t ); 
    y = 0.5 * height + R * sin ( t ); label = C4; };

int n = 7;

mesh Th = buildmesh
(
    floor   ( 5 * ( width / height ) * n ) 
  + right   ( 5 * n )
  + ceiling ( 5 * ( width / height ) * n )
  + left    ( 5 * n )
  + cir     ( 5 * n )
);

fespace Xh ( Th, P2 );
fespace Mh ( Th, P1 );
Xh u2, v2, u1, v1, up1, up2;
Mh p, q;
Xh v, u, w;
Xh uold = 0.0;
Xh vcty;

plot ( Th, wait = false,
  ps = "karman_vortex_mesh.ps" );

int nRuns = 700;
bool reuseMatrix = false;
real dt = 7.0 * scale * scale;
//
//  some fantasy heat conducting water to create a nice looking temperature plot
//
real vinf = 0.0018 / scale;  // m/s
real rho = 1000.0;           // density kg/m^3
real nu = 0.9e-6;            // kinematic viscosity m^2/s
real cp = 4.2 * 1000.0;      // heat capacity J/kg*K
real lambda = 7 * 0.6;       // heat conductivity W/mK
real mu = nu * rho;
real a = lambda / ( rho * cp );

cout << endl;
cout << "  Reynolds Number: " << 2.0 * R * vinf / nu << endl;
cout << "  Prandtl Number: " << nu / a << endl;

problem NS ( [u1,u2,p], [v1,v2,q], solver=UMFPACK, init=reuseMatrix ) =
  int2d(Th)
  ( 
    u1 * v1 + u2 * v2                      //  From advective term in Lagrange coordinates
    + dt*nu*(dx(u1)*dx(v1) + dy(u1)*dy(v1) //  Viscous heat
    + dx(u2)*dx(v2) + dy(u2)*dy(v2))
    + p*q*1.e-6                            //  Stabilization
    - dt*(p*(dx(v1) + dy(v2))/rho          //  Pressure term
    + q*(dx(u1) + dy(u2)))                 //  Mass conservation ( Variational continuity equation )
  )               
  - int2d ( Th )
  (
      convect ( [up1,up2], -dt, up1 ) * v1
    + convect ( [up1,up2], -dt, up2 ) * v2
  )
  + on ( C1, u2 = 0.0 )
  + on ( C2, u1 = vinf, u2 = 0.0 )
  + on ( C4, u1 = 0.0, u2 = 0.0 )
  + on ( C6, u2 = 0.0 )
;
//
//  Heat transfer.
//
real Tsurf = 1.0;

problem convectdiffusion ( u, v ) =
    int2d ( Th ) 
    (
      u * v + a * dt * ( dx ( u ) * dx ( v ) + dy ( u ) * dy ( v ) )
    )
  - int2d ( Th )
    (
      convect ( [ up1, up2 ], -dt, uold ) * v
    )
  + on ( C2, u = 0.0 )
  + on ( C4, u = Tsurf );
//
//  Write out data for later processing by MATLAB/Octave.
//
savemesh ( Th, "karman_vortex.msh" );
ffSaveVh ( Th, Mh, "karman_vortex_mh.txt" );
ffSaveVh ( Th, Xh, "karman_vortex_xh.txt" );

int k;
real t = 0.0;

for ( int i = 0; i <= nRuns ; i++ )
{
  up1 = u1;
  up2 = u2;
  NS;
  reuseMatrix = true;
  convectdiffusion;
  uold = u;
  k = i + 10000;
//
//  Temperature plot - clip to make the plot looking better
//
  w = u * ( abs ( u ) < 0.3 );

  plot ( w, wait = false, value = false, fill = true, ShowAxes = false, 
    cmm = "RunNumber: "+i+"/"+nRuns+"  Time: "+t+"sec"+"  vInf: "+vinf+"m/s" );
  t = t + dt;
//
//  Compute the vorticity.
//
  vcty = dx ( u2 ) - dy ( u1 );
//
//  Write time-dependent data to files for later processing by MATLAB/Octave.
//
  ffSaveData ( u, "karman_vortex_temperature_"+k+".txt" );
  ffSaveData ( vcty, "karman_vortex_vorticity_"+k+".txt" );
  ffSaveData ( p, "karman_vortex_pressure_"+k+".txt" );
}

  plot ( w, wait = false, value = false, fill = true, ShowAxes = false, 
    cmm = "RunNumber: "+i+"/"+nRuns+"  Time: "+t+"sec"+"  vInf: "+vinf+"m/s",
    ps = "karman_vortex.ps");
//
//  Terminate.
//
cout << "\n";
cout << "karman_vortex:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );