// cylinder.edp
//
//  Discussion:
//
//    Flow around a cylinder: fixed point iteration solver
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/examples/freefem++/cylinder.edp
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    15 December 2014
//
//  Author:
//
//    Hoang Tran
//

//string to save file-outputs
string fnse = "cylinder";

//Define the domain 
//rectangle [L,R] X [B,T]
real L=0.0,R=2.2;
real B=0.0,T=0.41;
int  nBottom1=45, nTop1=45,nBottom2=75,nTop2 = 75;
int  nRight=45,nLeft=30;
int  nc=180;

// outer rectangle
border Bottom1(t=0,0.50/R){x=L+(R-L)*t;y=B;label=1;};
border Top1(t=0,0.50/R){x=L+(R-L)*t;y=T;label=3;};
border Bottom2(t=0.50/R,1){x=L+(R-L)*t;y=B;label=11;};
border Top2(t=0.50/R,1){x=L+(R-L)*t;y=T;label=31;};
border Left(t=0,1){x=L;y=B+(T-B)*t;label=4;};
border Right(t=0,1){x=R;y=B+(T-B)*t;label=2;};
// obstacle
border Obs(t=0,2*pi){x=0.05*cos(t)+0.2; y=0.05*sin(t)+0.2;label=5;};

// build the mesh 
mesh ThF=buildmesh(Bottom1(nBottom1)+ Bottom2(nBottom2)+Right(nRight)+Top1(-nTop1) + Top2(-nTop2)+Left(-nLeft)
              +Obs(-nc));

// plot the mesh           
//plot(ThF,wait=1);
 
real dT = 0.005;
real T0 = 0.0;
real Tf = 8.0;         
real PRESSUREPENALTY = 1.0e-12;   
real ERRORTOL = PRESSUREPENALTY*1.0e2; 

//viscosity
real nu  = 1.0e-3;
//mean inflow velocity
real Umax= 1.0; 
//diameter of cylinder 
real diameter  = 0.1;
//density
real rho = 1.0;

int  TIMELEVELFREQUENCY = 200;   // number of steps between data-output
int  TOTALTIMELEVELS    = 8;     // (Tf - T0) / (dT * TIMELEVELFREQUENCY) 
int  MAXITS = 300;

//********************************************************************
//********************************************************************
//********************************************************************
//FE-spaces
//Velocity space, C0- piecewise quadratic
  fespace Uh(ThF,P2);
//  plot(ThF);
  
//Pressure space, C0- piecewise linear
  fespace Ph(ThF,P1);
  
//fe-functions (velocity space)
// uXold,  uYold  = previous time-step
// uXtemp, uYtemp = current time-step, previous fixed point iteration
// uX,     uY     = current time-step, current iteration
  Uh uX, uY, uXold, uYold, uXTemp, uYTemp;
  
//fe-functions (velocity space, test functions)
  Uh vX, vY;
  
  Uh kX, kY; 
  
//fe-functions (pressure space)
// p     = pressure at current time-step
// pold  = pressure at previous time-step
  Ph p, pold;
  
//fe-functions (pressure space, test functions)
  Ph q;
  
func real uXb(real x, real y, real t){
  return (6.0/(0.41)^2 *sin(pi*t/8.0)*y*(0.41-y));
  }
  
//Define expressions via macros:

macro div(u1,u2) ( dx(u1) + dy(u2) ) //
macro dot(u1,u2,v1,v2) ( u1 * v1 + u2 * v2 ) //
macro Ugradv1(u1,u2,v1) ( u1 * dx(v1) + u2 * dy(v1) ) //

macro cc(u1,u2,v1,v2,w1,w2) ( Ugradv1(u1,u2,v1) * w1 + Ugradv1(u1,u2,v2) * w2 ) //
macro cch(u1,u2,v1,v2,w1,w2) ( 0.5*(cc(u1,u2,v1,v2,w1,w2) - cc(u1,u2,w1,w2,v1,v2)) ) //  
macro contract(u1,u2,v1,v2) ( dx(u1) * dx(v1) + dx(u2) * dx(v2) + dy(u1) * dy(v1) + dy(u2) * dy(v2) ) //

//---------------------------------------------------------------------   

//********************************************************************
//********************************************************************
//********************************************************************

//---------------------------------------------------------------------    

//solver option
// plotsDEBUG = 0  --  do not produce plots of approximations
// plotsDEBUG = 1  --  print plots to screen
// plotsDEBUG = 2  --  output to data file (uX, uY, p, phi)
int   plotsDEBUG = 2;

//MODEL PARAMETERS 

//stochastic viscosity

//---------------------------------

// initial codition 

func real uXInit(real x, real y){
	return (0.0); 
}

func real uYInit(real x, real y){
	return (0.0); 
}

func real pInit(real x, real y){
	return (0.0);
}

//-------------------
//  Output flow details

{
    ofstream report(fnse + "_REPORT.txt");
    report << " ***            BEGIN REPORT, Hoang Tran              ***" << endl;
    report << " ***            Flow around a cylinder                ***" << endl << endl;
    report << "Numerical scheme: " << endl;
    report << " *** Crank-Nicolson scheme     ***" << endl;
    report << " *** Fixed point iteration     ***" << endl << endl;
    report << "PROBLEM SETUP: From Section 2, John, Int J Numer Meth Fluids, 2004" << endl;
    report << " # bdry nodes (channel length bottom), N01  = " << nBottom1 + nBottom2 << endl;  
    report << " # bdry nodes (channel height right),  N02  = " << nRight << endl;  
    report << " # bdry nodes (channel length top),    N03  = " << nTop1 + nTop2 << endl;  
    report << " # bdry nodes (channel height left),   N04  = " << nLeft << endl; 
    report << " # bdry nodes (obstacle),              N05  = " << nc << endl; 
    report << "TIME STEPPING DATA: " << endl;
    report << "  Time step size,               dT = " << dT << endl;
    report << "  Initial time,                 T0 = " << T0 << endl;
    report << "  Final time,                   Tf = " << Tf << endl;
    report << "PROBLEM PARAMETERS: " << endl;
    report << "Fluid viscosity,         nu        = " << nu << endl << endl;
}

//*********** Main SOLVE (crank-nicolson, fixed point iterations) ****************
	
//****************************************************************************************************************************************
//****************************************************************************************************************************************
//****************************************************************************************************************************************

//initialize
// temp1, temp2     = temporary holder for error at each nonlinear iteration
//                initialized to enter loop
real temp1, temp2;
 
// initialize 
// time t at T0
// time index tracker timeIndex 
//   and timeString to store timeIndex as a string
real   t = T0;
int    timeIndex = 1;
string timeString;

problem cnITER([uX,uY,p],[vX,vY,q],solver=UMFPACK) =
   int2d(ThF)(
      dot(uX,uY,vX,vY)
    + 0.5 * dT * cch( uXTemp, uYTemp, uX, uY, vX, vY)
    + 0.5 * dT * nu * contract(uX, uY, vX, vY)
    - dT * q * div(uX, uY)
    - dT * p * div(vX, vY)
    + dT * PRESSUREPENALTY * p * q
   )
 + int2d(ThF)( 
    - dot(uXold,uYold,vX,vY)   
    + 0.5 * dT * cch( uXold, uYold, uXold, uYold, vX, vY)
    + 0.5 * dT * nu * contract(uXold, uYold, vX, vY)
   )
 +on( 2,4, uX = uXb(x,y,t) , uY = 0.0 )
 +on( 1,11,3,31,5, uX = 0.0 , uY = 0.0 )  ;
 
//Get test function Vd, Vl for drag and lift calculation
Uh phi1, phi2, phiTest1, phiTest2;  //Vd = (phi1, phi2)  Vl = (phi2, phi1)
func g = 0.0;
problem getPhi ([phi1,phi2], [phiTest1,phiTest2], solver = UMFPACK) =  //get phi1 such that phi1 = 1 on 5 & 0 on others
	int2d(ThF)(
		contract(phi1,phi2,phiTest1,phiTest2)
	)
	- int2d(ThF)( //this is probably unnecessary
		g*phiTest1 + g*phiTest2
	)
	+ on(5, phi1 = 1.0, phi2 = 0.0)  //phi2 = 0 on all borders
	+ on(1,2,3,4,11,31, phi1 = 0.0, phi2 = 0.0);
getPhi;

real kinetic, dissipation, dragCoef, liftCoef;

real relError;
int  numIts, numTotalIts;

// initialize
// Take (u(0),p(0)) (initial guess) to be the solution to the stokes problem
// OUTPUT:
//   uX, uY, p
uX  = uXInit(x,y);      
uY  = uYInit(x,y);
p   = pInit(x,y);
    
//initialize real-time clock tracker
real outputCLOCK =   clock();

//LOOP #1:  Time evolution for nse
while ( t < (Tf-1.0e-12) ){

  cout << endl << "NUMERICAL TIME LEVEL  t = " << t << endl;	
  cout << endl << " *** real time    = " << clock() << endl;	
  cout << endl << " *** elapsed time = " << clock() - outputCLOCK << endl;
  outputCLOCK = clock();
  
  //Update time-level
  t = t + dT;

  //Store from previous time-step, 
  //  u(n) <- u(n+1)
  uXold  = uX;      
  uYold  = uY;
  pold   = p;
  
  //initialize
  // temp     = temporary holder for error at each nonlinear iteration
  //                initialized to enter loop
  // relError = temporary holder for relative error at each nonlinear iteration
  //            -- initialized to enter loop
  // numIts   = iteration tracker for nonlinear loop
  relError  = 1.0;
  numIts    = 0;
  
  //LOOP #1a: Fixed point iterations for nonlinear solve
  while ( (relError > ERRORTOL) & (numIts < MAXITS)) {
   
  //update iteration number
  numIts++;
  numTotalIts++;

   //Store from previous newton iterate, 
   //  u(n+1,j) <- u(n+1,j+1)
   uXTemp  = uX;      
   uYTemp  = uY;
    
   //Solve nse at current time level
   //OUTPUT:
   //   uX(n+1,j+1), uY(n+1,j+1), p(n+1,j+1)
   cnITER;
   
   kX = uXTemp - uX;
   kY = uYTemp - uY;
         
   //Compute error in iteration
   temp1 = sqrt( int2d(ThF)(  (kX)^2     + (kY)^2 
                            + (dx(kX))^2 + (dy(kY))^2 
                            + (dy(kX))^2 + (dx(kY))^2 
                            )); 
                            
   temp2 = sqrt( int2d(ThF)(  (uXTemp)^2     + (uYTemp)^2 
                            + (dx(uXTemp))^2 + (dy(uYTemp))^2 
                            + (dy(uXTemp))^2 + (dx(uYTemp))^2 
                            )); 
                            
    if (numIts == 1){
      relError = 1.0;
    }
    else{
      relError = temp1 / temp2;
    }
                                
  } // end of while loop #1a (newton/picard iteration)
  
  if ( relError > ERRORTOL ){
   {
    ofstream report(fnse + "_REPORT.txt",append);
    report << endl << " ***********************************************************************  " << endl;
    report << " **** WARNING: nonlinear iterates DO NOT CONVERGE,  t = " << t << endl << endl;
    report << "      numIts        = " << numIts << endl;
    report << "      error(u)    = " << temp1 << endl;
    report << "      relError(u) = " << relError << endl;
    report << endl << " ***********************************************************************  " << endl;
   } 
  }      
  
  cout << endl << "Nonlinear solve, iteration:  time = " << t << endl;
  cout << "  numIts       = " << numIts << endl;
  cout << "  error(u)   = " << temp1 << endl;
  cout << "  relError(u)= " << relError << endl;
  
  if( timeIndex == TIMELEVELFREQUENCY ){
    timeString = t;
    timeIndex = 1;
    
    if ( plotsDEBUG == 1 ){
      plot(p, [uX,uY], fill=1, cmm="Velocity, pressure fields:  time level = " + timeString,wait=1,value=1);
    }
    else if (plotsDEBUG == 2) {
//	  {
//      plot(p,[uX,uY], value=1, fill=1, ps=fnse + "_" + timeString + ".eps"); 
//      } 
      {
        ofstream file(fnse + "_velocityX_TIME_" + timeString + ".txt");
        file << uX[] << endl;
      }
      {
        ofstream file(fnse + "_velocityY_TIME_" + timeString + ".txt");
        file << uY[] << endl;
      }
      {
        ofstream file(fnse + "_pressure_TIME_" + timeString + ".txt");
        file << p[] << endl;
      }
    }
  }
  else{
     timeIndex++;
  }
  
//**********************************************************************
//**********************************************************************   
//COMPUTE PRESSURE DIFF & DRAG AND LIFT COEFFICIENTS 

  dragCoef = -20.0*int2d(ThF)( (1.0/dT)*dot(uX,uY,phi1,phi2) - (1.0/dT)*dot(uXold,uYold,phi1,phi2) + nu*contract(uX,uY,phi1,phi2) + cc(uX,uY,uX,uY,phi1,phi2) - p*div(phi1,phi2) );
  liftCoef = -20.0*int2d(ThF)( (1.0/dT)*dot(uX,uY,phi2,phi1) - (1.0/dT)*dot(uXold,uYold,phi2,phi1) + nu*contract(uX,uY,phi2,phi1) + cc(uX,uY,uX,uY,phi2,phi1) - p*div(phi2,phi1) );
  kinetic  =  0.5*sqrt( int2d(ThF)(  (uX)^2     + (uY)^2));
  dissipation = 0.5*sqrt( int2d(ThF)( (dx(uX))^2 + (dy(uY))^2 
                            + (dy(uX))^2 + (dx(uY))^2 ));    


 {
        ofstream file(fnse + "_DRAG.txt",append);
        file << t << "     " << dragCoef << endl;
  }

  {
        ofstream file(fnse + "_LIFT.txt",append);
        file << t << "     " << liftCoef << endl;
  }

  {
        ofstream file(fnse + "_KINETICS.txt",append);
        file << t << "     " << kinetic << endl;
  }

  {
        ofstream file(fnse + "_DISSIPA.txt",append);
        file << t << "     " << dissipation << endl;
  }
  
} // END WHILE LOOP #1 (time stepping)

{
   ofstream report(fnse + "_REPORT.txt",append);
   report << endl << endl;
   report << "***********************************************************************" << endl;
   report << "  Total # nonlinear-iterations     = " << numTotalIts << endl;
   report << "  Total time                       = " << clock() << endl;
 }