//######################################################################################
//
// Copyright, 2009, Wang Xian, Tokyo Institute of technology                                           
// -- HSMAC.c (Feb.2009)                                                                                         
// -- Using HSMAC (Highly Simplified Marker and Cell) method to solve the 2D cavity flow 
// -- CPU code
                                                                                                   
//######################################################################################
                                 
                                                                
     #include <stdio.h>
     #include <stdlib.h>
     #include <malloc.h>
     #include <assert.h> //what is this?
     #include <math.h>
     #include <time.h>
     #include <sys/time.h>
     #include <string.h>

     #define DOUBLE float 


     #ifndef max
     #define max(x,y) (((x)<=(y)) ? (y):(x))
     #endif

     #ifndef min
     #define min(x,y) (((x)>=(y)) ? (y):(x))
     #endif

     #define abs(x) (((x)<0) ? (-x):(x)) /* NOTICE: when using expression a+b, abs((a+b)) should be used instead of abs(a+b) */


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

void init (int,int,DOUBLE,DOUBLE**,DOUBLE**,DOUBLE**);

DOUBLE HSMAC_CPU(int,int,DOUBLE,DOUBLE,DOUBLE,DOUBLE,DOUBLE,DOUBLE,DOUBLE,DOUBLE**,DOUBLE**,DOUBLE**,DOUBLE**,DOUBLE**);

void final_output(int,int,int,DOUBLE,DOUBLE,DOUBLE**,DOUBLE**,DOUBLE**) ;

void computation_con(int,int,DOUBLE**,DOUBLE**,DOUBLE**);

DOUBLE cpu_time();



//********************** START MAIN ********************************************************************

int main(int argc, char *argv[])
{   
	
	  char  input[20];                       

	  int   nx, ny, cycle, cycle_max, output, tt, yy;   // real nodes: 1--(nx-2)  1--(ny-2)
                                                        // imaginary nodes: (0 & nx-1) (0 & ny-1)	  
	  DOUBLE dt,dx,dy,lx,ly;

	  DOUBLE Re,u_0,omg,epsp;
  
	  DOUBLE **p0, **u0, **v0, **un, **vn; 

	  DOUBLE  cpu_t=0.0, cpu_t0=0.0, delt_t=0.0, divmax_cpu=0.0;
    	  
      FILE  * time_data, * in_data;  
	

//-------------------------------------------------------------------------	
	if((in_data=fopen("INPUT.dat","r"))==NULL)
	{
		printf("cannot open this file\n");
		exit(0);
	}


	if((time_data=fopen("time.dat","w"))==NULL)
	{
		printf("cannot open this file\n");
		exit(0);
	}


//read paramenter:	nx,ny,cycle,cycle_max,output,Re,rho_0,u_0;

	 fscanf(in_data,"%s%d\n",input,&nx);
	 fscanf(in_data,"%s%d\n",input,&ny);
	 fscanf(in_data,"%s%d\n",input,&cycle);                /* time step, that is loop counter */
	 fscanf(in_data,"%s%d\n",input,&cycle_max);            /* maximum number of loop */ 
	 fscanf(in_data,"%s%d\n",input,&output);              /* output interval */

	 /* NOTICE: double type parameter should use %lf or %le to read !!*/
	 fscanf(in_data,"%s%f\n",input,&Re);                
	 fscanf(in_data,"%s%f\n",input,&dt);                
	 fscanf(in_data,"%s%f\n",input,&omg);      
	 fscanf(in_data,"%s%f\n",input,&epsp); 
	
	 fclose(in_data);
	 
    printf("-----Please check the input data!!-----\n");
    printf("nx=%d ny=%d cycle=%d cyclemax=%d output=%d\n",nx,ny,cycle,cycle_max,output);
    printf("Re=%f dt=%f omg=%f epsp=%f\n",Re,dt,omg,epsp);

//-------------------------------------------------------------------------------
	lx = 1.0;
	ly = 1.0;
    u_0 = 1.0; //(non-dimensional)
	dx = lx/(float)(nx-2);
	dy = ly/(float)(ny-2);


//----------------------------------------------------------------------------------
    u0  = NULL;
    v0  = NULL;
    p0  = NULL;
    un  = NULL;
    vn  = NULL;

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

	/* allocate 2D density rho */
	u0=(DOUBLE **)calloc(ny,sizeof(DOUBLE *));
    for(yy=0;yy<ny;yy++)u0[yy]=(DOUBLE *)calloc(nx,sizeof(DOUBLE));
    
	v0=(DOUBLE **)calloc(ny,sizeof(DOUBLE *));
    for(yy=0;yy<ny;yy++)v0[yy]=(DOUBLE *)calloc(nx,sizeof(DOUBLE));

	p0=(DOUBLE **)calloc(ny,sizeof(DOUBLE *));
    for(yy=0;yy<ny;yy++)p0[yy]=(DOUBLE *)calloc(nx,sizeof(DOUBLE));

	un=(DOUBLE **)calloc(ny,sizeof(DOUBLE *));
    for(yy=0;yy<ny;yy++)un[yy]=(DOUBLE *)calloc(nx,sizeof(DOUBLE));

	vn=(DOUBLE **)calloc(ny,sizeof(DOUBLE *));
    for(yy=0;yy<ny;yy++)vn[yy]=(DOUBLE *)calloc(nx,sizeof(DOUBLE));

//---------------------------------------------------------------------------------	
	if(cycle==0)
	{
      init(nx,ny,u_0,p0,un,vn);
      tt=1;
  	}
	else
	{
      computation_con(nx,ny,p0,un,vn);
      tt=cycle;
	}

//######################### START LOOP #############################################	

	 cpu_t0 = cpu_time();                  /* start to record the time */
	 cpu_t = cpu_t0;

	 {
        for(cycle=tt;cycle<=cycle_max;cycle++)
		{
          divmax_cpu=HSMAC_CPU(nx,ny,u_0,dx,dy,dt,omg,epsp,Re,u0,v0,p0,un,vn);
        
			if((cycle%output)==0) 
		    {
			  delt_t = cpu_time()-cpu_t;      /* time for every "output" time steps */
	          cpu_t = cpu_time()-cpu_t0;      /* total time until now */

			  printf("----timestep = %d   divmax=%g---------------\n",cycle,divmax_cpu);

              final_output(nx,ny,cycle,dx,dy,p0,un,vn);
	          
			   fprintf(time_data,"cycle=%d\n",cycle);
			   fprintf(time_data,"divmax=%g\n",divmax_cpu);
			   fprintf(time_data,"cputime=%g, delta_t=%g\n",cpu_t,delt_t);
	           fprintf(time_data,"---------------------------------------------\n");
			}
		}
	 }
	 

		 		  
//########################### END OF LOOP ################################################

	fclose(time_data);

	return 0;
}/***** end main *****/



//************** init() ********************************
void init (int nx,int ny,DOUBLE u_0,DOUBLE **p0,DOUBLE **un,DOUBLE **vn)
{
   int x,y;
	   
	   for(y=0;y<ny;y++)
	   {
		   for(x=0;x<nx;x++)
		   {
			   p0 [y][x] = 0.0;                 
			   un [y][x] = 0.0;
			   vn [y][x] = 0.0;
	           
			   if(y==ny-1) un[y][x]=2.0*u_0-un[y-1][x];          /* the top is moving lid with u_0 */

		    }
	    }


}/***** end init *****/


//****************** HSMAC_CPU() ********************************

DOUBLE HSMAC_CPU(int nx,int ny,DOUBLE u_0,DOUBLE dx,DOUBLE dy,DOUBLE dt,DOUBLE omg,DOUBLE epsp,DOUBLE Re,
	             DOUBLE **u0,DOUBLE **v0,DOUBLE **p0,DOUBLE **un,DOUBLE **vn)
{   
  int x,y,flag,iter;
  DOUBLE uu,vv,fux,fuy,difu,fvx,fvy,difv;
  DOUBLE divmax1 = 0.0, delp, del, div;
		
  for(y=0;y<ny;y++)
  {
    for(x=0;x<nx;x++)
	{

	   u0 [y][x] = un [y][x] ; 
	   v0 [y][x] = vn [y][x] ;

	}
  }


  //------------------ velocity field ---------------------------------------------------------
  
  
  for(y=1;y<ny-1;y++)
  {
    for(x=1;x<nx-1;x++)
	{
   
//---------------------- compute u ----------------

	   vv = (v0[y][x]+v0[y][x+1]+v0[y-1][x+1]+v0[y-1][x]) / 4.0 ;
       
	   fux = u0[y][x] * (u0[y][x+1]-u0[y][x-1])/(2.*dx) ;
	 
       fuy =    vv * (u0[y+1][x]-u0[y-1][x])/(2.*dy) ;  
  
       difu= 1./Re * (  (u0[y][x-1]-2.0*u0[y][x]+u0[y][x+1])/(dx*dx) 
                       +(u0[y-1][x]-2.0*u0[y][x]+u0[y+1][x])/(dy*dy)  ) ;

	   un[y][x] = u0[y][x] + dt * (-fux - fuy + difu + (p0[y][x]-p0[y][x+1])/dx ) ;
   
//---------------------- compute v ----------------

	   uu = (u0[y][x]+u0[y][x-1]+u0[y+1][x-1]+u0[y+1][x]) / 4.0 ;
       
	   fvx =    uu * (v0[y][x+1]-v0[y][x-1])/(2.*dx) ;
	 
       fvy = v0[y][x] * (v0[y+1][x]-v0[y-1][x])/(2.*dy) ;  
  
       difv= 1./Re * (  (v0[y][x-1]-2.0*v0[y][x]+v0[y][x+1])/(dx*dx) 
                       +(v0[y-1][x]-2.0*v0[y][x]+v0[y+1][x])/(dy*dy)  ) ;

	   vn[y][x] = v0[y][x] + dt * (-fvx - fvy + difv + (p0[y][x]-p0[y+1][x])/dy ) ;

	}
  }



//*----------------- velocity BC --------------------

  for(y=0;y<ny;y++)
  {  
	  un[y][0] = 0.0 ; 
	  un[y][nx-2] = 0.0 ; 

	  vn[y][0] = -vn[y][1] ;
      vn[y][nx-1] = -vn[y][nx-2] ;
  }


  for(x=0;x<nx;x++)
  {  
	  un[0][x] = -un[1][x] ;
	  un[ny-1][x] = 2.0*u_0-un[ny-2][x] ;

	  vn[0][x] = 0.0 ;
      vn[ny-2][x] = 0.0 ;
  }

  
  //--------------- pressure field ---------------------------*/

flag=0;
iter=0;

loop: divmax1=0;
for(y=1;y<ny-1;y++)
  {
    for(x=1;x<nx-1;x++)
	{

	   del = dt * (2.0/(dx*dx) + 2.0/(dy*dy) ) ;  
       
	   div = ( un[y][x]-un[y][x-1] )/dx + ( vn[y][x]-vn[y-1][x] )/dy ;
  
	   divmax1 = max(divmax1,abs(div)) ;

	   delp = -omg * div / del ;

	   p0[y][x] = p0[y][x] + delp ;
	   
	   un[y][x] = un[y][x] + dt/dx * delp ;
	   
	   un[y][x-1] = un[y][x-1] - dt/dx * delp ;

	   vn[y][x] = vn[y][x] + dt/dy * delp ;
	   
	   vn[y-1][x] = vn[y-1][x] - dt/dy * delp ;
	}
  }
	    
	    
  for(y=0;y<ny;y++)
  {  
	  un[y][0] = 0.0 ; 
	  un[y][nx-2] = 0.0 ; 

	  vn[y][0] = -vn[y][1] ;
      vn[y][nx-1] = -vn[y][nx-2] ;
  }


  for(x=0;x<nx;x++)
  {  
	  un[0][x] = -un[1][x] ;
	  un[ny-1][x] = 2.0*u_0-un[ny-2][x] ;

	  vn[0][x] = 0.0 ;
      vn[ny-2][x] = 0.0 ;
  }
	    
	    
	    
	    iter = iter +1;

	   if ( (divmax1 < epsp) || iter == 1000) flag=1;

	   if (flag == 0) goto loop;


return(divmax1);

}/***** end of HSMAC_CPU *****/




//*********************************** final_output() ************************************
void final_output(int nx,int ny,int cycle,DOUBLE dx,DOUBLE dy,DOUBLE **p0,DOUBLE **un,DOUBLE **vn) 
{

	int x,y;
	int kk1,kk2,kk3,kk4,kk5,kk6,kk7;
	int kk8,kk9,kk10,kk11,kk12,kk13,kk14,kk15,kk16;
	DOUBLE xx,yy,uu,vv,pp;


	FILE  * file_dat;
	FILE  * copy_dat;
//	FILE  * fp;

	char data_name[200];


    kk1=cycle;
	kk2=kk1/10000000;
	kk3=kk1-kk2*10000000;
	kk4=kk3/1000000;
	kk5=kk3-kk4*1000000;
	kk6=kk5/100000;
	kk7=kk5-kk6*100000;
	kk8=kk7/10000;
	kk9=kk7-kk8*10000;
	kk10=kk9/1000;
	kk11=kk9-kk10*1000;
	kk12=kk11/100;
	kk13=kk11-kk12*100;
	kk14=kk13/10;
	kk15=kk13-kk14*10;
	kk16=kk15/1;



//------------------------------------- non-binary tecplot data file-----------------------------------

	sprintf(data_name,"%d%d%d%d%d%d%d%d.plt",kk2,kk4,kk6,kk8,kk10,kk12,kk14,kk16);
	
	file_dat=fopen(data_name,"w");

	fputs("VARIABLES=X,Y,U,V,P\n",file_dat);
	fprintf(file_dat, "ZONE I=%6d, J=%6d, F=POINT\n",nx-1,ny-1);

		for (y=0; y<ny-1; y++) 
		for (x=0; x<nx-1; x++) 
		{{
			xx=dx*x;
			yy=dy*y;
			uu=(un[y][x]+un[y+1][x])/2.0;
			vv=(vn[y][x]+vn[y][x+1])/2.0;
			pp=(p0[y][x]+p0[y][x+1]+p0[y+1][x]+p0[y+1][x+1])/4.0;
		   fprintf(file_dat,"%e %e %20.7e %20.7e %20.7e\n",xx,yy,uu,vv,pp);
	    }
	}

	fclose(file_dat);

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


// copy data for continue the computation in the future

	if((copy_dat=fopen("continue.dat","w"))==NULL)
	{
		printf("cannot open this file\n");
		exit(0);
	}


	for(y=0; y<ny; y++)
	{
	    for(x=0; x<nx; x++)
	    {
		   fprintf(copy_dat,"%e %e %e\n",p0[y][x],un[y][x],vn[y][x]); 
	    }
	}

	fclose(copy_dat);


}/***** end final_output *****/



//*********************************** computation_con() ************************************
void computation_con(int nx,int ny,DOUBLE **p0,DOUBLE **un,DOUBLE **vn)
{
	int x,y;
	FILE  * read_dat;


	if((read_dat=fopen("continue.dat","r"))==NULL)
	{
		printf("cannot open this file\n");
		exit(0);
	}


	for(y=0; y<ny; y++)
	{
	    for(x=0; x<nx; x++)
	    {
  	       /* NOTICE: double type parameter should use %lf or %le to read !!*/
		   
	    fscanf(read_dat,"%e %e %e\n",&p0[y][x],&un[y][x],&vn[y][x]);
	    }
	}
	
	fclose(read_dat);

}/***** end computation_con *****/


//*************************************** cpu_time()**********************************
DOUBLE cpu_time()
{

	 struct timeval tm;
     DOUBLE t ;
     static int base_sec = 0,  base_usec = 0;

     gettimeofday(&tm, NULL);
  
     if(base_sec == 0 && base_usec == 0) 
	 {
        base_sec = tm.tv_sec;   
		base_usec = tm.tv_usec;
        t = 0.0;
     }
     else 
	 {
        t = (DOUBLE) (tm.tv_sec-base_sec) + 
            ((DOUBLE) (tm.tv_usec-base_usec))/1.0e6 ;
     }

     return t ;

}/***** end cpu_time *****/