static char help[]="Solidification of a alloy \nAuther:Sepideh Kavousi PHD Student\nStart date: April 2 -2018\nIt is Model in Echebarria, floch paper for directional solidification, isotropic \n in this file I did not consider antitrapping (U2=0) and try to converge results";

#include <petsc.h>
#include "common.h"

// one_%f-%06d.vts: pressure, MAC, ibrun -n 68 ./one.out -ts_monitor -snes_fd_color -ts_max_snes_failures -1 -ts_type bdf -ts_bdf_adapt -pc_type bjacobi -snes_linesearch_type l2 -snes_type ksponly -ksp_type gmres -ksp_gmres_restart 1001 -sub_pc_type bjacobi -sub_ksp_type preonly
// this did not work






//old written 
// mpiexec -n 4 ./one.out  -ts_monitor -snes_monitor -snes_fd_color -ts_max_snes_failures -1  -ts_type arkimex -ts_adapt_type none -snes_type newtontr -pc_type bjacobi -snes_linesearch_monitor -snes_linesearch_type bt -snes_type newtontr -ksp_type gmres

//options for run on mike
///work/skavou1/petsc/one/bin/mpiexec -n 16 ./one.out -ts_monitor -snes_monitor -snes_fd_color -ts_max_snes_failures -1  -ts_type arkimex -ts_adapt_type none -snes_type newtontr -pc_type bjacobi -snes_linesearch_monitor -snes_linesearch_type bt -snes_type newtontr -ksp_type gmres

//options from Moose
//./one.out -ts_monitor -snes_fd_color -ts_max_snes_failures -1  -ts_type arkimex -ts_adapt_type none -snes_type newton -pc_type asm  -snes_linesearch_type bt -snes_type newtontr -ksp_type gmres -ksp_gmres_restart 31 -sub_ksp_type preonly -sub_pc_type lu -pc_asm_overlap 2


//options for fastest
//./one.out -ts_monitor -snes_fd_colar -ts_max_snes_failures -1  -ts_type beuler -pc_type bjacobi  -snes_linesearch_type nleqerr -snes_type newtontr -ksp_gmres_restart 31 -sub_pc_type ilu -snes_monitor -snes_converged_reason -ksp_monitor


//options for anisotropic
//./one.out -ts_monitor -snes_fd_colar -ts_max_snes_failures -1  -ts_type beuler -pc_type bjacobi  -snes_linesearch_type bt -snes_type newtontr -ksp_gmres_restart 31 -sub_pc_type ilu -snes_monitor -snes_converged_reason -ksp_monit -ksp_type gmres

typedef struct {
double U,p,vy,vx,pp;
} Field;

//initial boundary condition
PetscErrorCode InitialState(Vec y,DMDALocalInfo *info2,struct VAR_STRUCT *user) {
	PetscErrorCode ierr;
	int j,i;
	Field **aY;

	ierr= DMDAVecGetArray(user->da,y,&aY) ;CHKERRQ(ierr);
ierr= PetscPrintf(PETSC_COMM_WORLD,"initial!\n") ;CHKERRQ(ierr);
	 for(j=info2->ys;j<info2->ys+info2->ym;j++){
		for(i=info2->xs;i<info2->xs+info2->xm;i++){
			aY[j][i].p=user->p_l;
			aY[j][i].U=user->U_l;
			aY[j][i].vy=user->vy_b;
			aY[j][i].vx=user->vx_b;
			aY[j][i].pp=0.0;
			if ((j>=user->Ny_Ib) & (j<=user->Ny_It)){
				if ((i>=user->Nx_Il) & (i<=user->Nx_Ir)){
					aY[j][i].p=user->p_s;
					aY[j][i].U=user->U_s;
					aY[j][i].vy=0.0;
					aY[j][i].vx=0.0;
					aY[j][i].pp=0.0;
					}

				}

			if ((i==0) && (j !=0) && (j !=info2->my-1)){
				aY[j][i].vy=user->vy_b;
				aY[j][i].vx=user->vx_b;
				aY[j][i].pp=user->pp_b;}



			}
		}

	ierr= DMDAVecRestoreArray(user->da,y,&aY) ;CHKERRQ(ierr);
	return 0;
	}









//Function Calculation
PetscErrorCode FormFunction(DMDALocalInfo *info2,PetscReal t,Field **aY,Field **aYdot,Field **aF, struct VAR_STRUCT *user) {
	PetscErrorCode ierr;
	int i,j,j_m,size;
	PetscScalar hx2,hy2,hx3,hy3,hxy,hx2y,hxy2,Px,Py,Pxx,Pyy,Pxy,Pxt,Pyt,Ux,Uxx,Uy,Uyy;
	PetscScalar Px_old,Py_old,Pxx_old,Pyy_old,Pxy_old,Pxxx_old,Pyyy_old,Pxxy_old,Pxyy_old,Uy_old;
	PetscScalar VVx_x,VVx_y,VVx_xx,VVx_yy,VVy_x,VVy_y,VVy_xx,VVy_yy,grad_VVx,grad_VVy,pp_x,pp_y,pp_xx,pp_yy;
	PetscScalar theta,thetax,thetay,AA;
	int RR,R,L,LL,B,BB,T,TT;
	PetscScalar e,ep,epp;
	PetscScalar U1,U2,num,num2;
	PetscScalar P1,grad_2,grad;
	PetscScalar a1,a2,a3,a4,a5,a6,a7,a8,a9;
	PetscScalar x1,x2,x3,x4,x5,x6,x7,x8,x9;
	PetscScalar y1,y2,y3,y4,y5,y6,y7,y8,y9;
	Field **asol_old;
	PetscScalar **arand;
	PetscReal st;


	TSGetTime(user->ts,&t);  
//			ierr= PetscPrintf(PETSC_COMM_WORLD,"start time= %f!\n",t) ;CHKERRQ(ierr);
	if (t>user->time){  
		user->time=t;
		ierr= PetscPrintf(PETSC_COMM_WORLD,"copy!\n") ;CHKERRQ(ierr);
		ierr= VecCopy(user->sol,user->sol_old); CHKERRQ(ierr);
		ierr= VecSetRandom(user->rand,user->r);CHKERRQ(ierr);		
		ierr= TSGetTimeStep(user->ts,&st);CHKERRQ(ierr);
		user->tw+=st;


		if (user->t_write % user->t_write_freq ==0) {
			ierr= PetscPrintf(PETSC_COMM_WORLD,"Write output at step= %d!\n",user->t_write) ;CHKERRQ(ierr);

			ierr= PetscFPrintf(PETSC_COMM_WORLD,user->fp," %d\t%f\n",user->t_write,t);

			ierr= PetscSNPrintf(user->filename,sizeof(user->filename),"fieldsplit_original_%f-%06d.vts",user->accelerate,user->t_write);CHKERRQ(ierr);
			ierr= TSMonitorSolutionVTK(user->ts,user->t_write,t,user->sol_old,&user->filename);CHKERRQ(ierr);}
		ierr= DMGlobalToLocalBegin(user->da,user->sol_old,INSERT_VALUES,user->sol_local);
		ierr= DMGlobalToLocalEnd(user->da,user->sol_old,INSERT_VALUES,user->sol_local);
		user->t_write+=1;
		}

	hx2=user->hx*user->hx;
	hy2=user->hy*user->hy;
	hxy=user->hx*user->hy;
	hx3=user->hx*user->hx*user->hx;
	hy3=user->hy*user->hy*user->hy;
	hx2y=user->hx*user->hx*user->hy;
	hxy2=user->hx*user->hy*user->hy;
//			ierr= PetscPrintf(PETSC_COMM_WORLD,"s1!\n") ;CHKERRQ(ierr);
	ierr= DMDAVecGetArray(user->da,user->rand,&arand) ;CHKERRQ(ierr);
	ierr= DMDAVecGetArray(user->da,user->sol_local,&asol_old) ;CHKERRQ(ierr);
//			ierr= PetscPrintf(PETSC_COMM_WORLD,"s2!\n") ;CHKERRQ(ierr);
	for (j=info2->ys;j<info2->ys+info2->ym;j++){
		for (i=info2->xs;i<info2->xs+info2->xm;i++){
//			ierr= PetscPrintf(PETSC_COMM_WORLD,"start time %f (%d,%d)!\t",t,i,j) ;CHKERRQ(ierr);

			//if (i==0) {i=0;}//aF[j][i].U=-aY[j][info2->mx-1].U+aY[j][i].U ; aF[j][i].p=-aY[j][info2->mx-1].p+aY[j][i].p ;aF[j][i].vy=-aY[j][info2->mx-1].vy+aY[j][i].vy ;aF[j][i].vx=-aY[j][info2->mx-1].vx+aY[j][i].vx;aF[j][i].pp=-aY[j][info2->mx-1].pp+aY[j][i].pp;}
			//else if (i==info2->mx-1) {i=info2->mx-1;}//aF[j][i].U=-aY[j][0].U+aY[j][i].U ; aF[j][i].p=-aY[j][0].p+aY[j][i].p ;aF[j][i].vy=-aY[j][0].vy+aY[j][i].vy ;aF[j][i].vx=-aY[j][0].vx+aY[j][i].vx;aF[j][i].pp=-aY[j][0].pp+aY[j][i].pp;}
			/*else*/ if (j==0) {aF[j][i].U=aY[j+1][i].U-aY[j][i].U; aF[j][i].p=aY[j+1][i].p-aY[j][i].p ;aF[j][i].vy=aY[j][i].vy;aF[j][i].vx=aY[j+1][i].vx-aY[j][i].vx;aF[j][i].pp=aY[j+1][i].pp-aY[j][i].pp;}//(aY[j][i].vy-user->vy_b)*0.5*(1-aY[j][i].p);}
			else if (j==info2->my-1) {aF[j][i].U=aY[j][i].U-user->U_l; aF[j][i].p=aY[j][i].p-user->p_l ;aF[j][i].vy=aY[j][i].vy;aF[j][i].vx=aY[j][i].vx;aF[j][i].pp=aY[j][i].pp;}//(aY[j][i].vy-user->vy_b)*0.5*(1-aY[j][i].p);}

			else {

	
				Px=(aY[j][i+1].p-aY[j][i-1].p)/(2.0*user->hx);
				Py=(aY[j+1][i].p-aY[j-1][i].p)/(2.0*user->hy);
				Pxx=((aY[j][i+1].p+aY[j][i-1].p-2.0*aY[j][i].p)/(hx2));
				Pxy=((aY[j+1][i+1].p+aY[j-1][i-1].p-aY[j+1][i-1].p-aY[j-1][i+1].p)/(4.0*hxy));
				Pyy=((aY[j+1][i].p+aY[j-1][i].p-2.0*aY[j][i].p)/hy2);
				Pxt=(aYdot[j][i+1].p-aYdot[j][i-1].p)/(2.0*user->hx);
				Pyt=(aYdot[j+1][i].p-aYdot[j-1][i].p)/(2.0*user->hy);
				if ((Px*Px)>0){
					theta=PetscAtanReal(Py/Px);
					}
				else {
					theta=PETSC_PI*0.5;
					}


				if (((Px*Px)+(Py*Py))>0.0 ) { //
					thetax=(Px*Pxy-Py*Pxx)/((Px*Px)+(Py*Py));
					thetay=(Px*Pyy-Py*Pxy)/((Px*Px)+(Py*Py));
					}
				else  { 
					thetax=0.0;
					thetay=0.0;
					}
				e=(1.0+user->epsilon*cos(4.0*theta));;
				ep=-(4.0*user->epsilon*sin(4.0*theta));
				epp=-(4.0*user->epsilon*4.0*cos(4.0*theta));
				grad=(4.0*(aY[j][i+1].U+aY[j][i-1].U+aY[j+1][i].U+aY[j-1][i].U)+1.0*(aY[j+1][i+1].U+aY[j+1][i-1].U+aY[j-1][i+1].U+aY[j-1][i-1].U)-20.0*aY[j][i].U)/(6.0*hx2);
				grad_2=(4.0*(aY[j][i+1].p+aY[j][i-1].p+aY[j+1][i].p+aY[j-1][i].p)+1.0*(aY[j+1][i+1].p+aY[j+1][i-1].p+aY[j-1][i+1].p+aY[j-1][i-1].p)-20.0*aY[j][i].p)/(6.0*hx2);

				Ux=(aY[j][i+1].U-aY[j][i-1].U)/(2.0*user->hx);
				Uy=(aY[j+1][i].U-aY[j-1][i].U)/(2.0*user->hy);
				Uxx=((aY[j][i+1].U+aY[j][i-1].U-2.0*aY[j][i].U)/(hx2));
				Uyy=((aY[j+1][i].U+aY[j-1][i].U-2.0*aY[j][i].U)/hy2);

				U1=user->D*user->tau_0/(user->W*user->W);
				U2= ((Px*Px+Py*Py)>0.0) ? user->a_p/sqrt(Px*Px+Py*Py) : 0.0 ;
				P1= (user->W*user->hy*(j-user->Ny_It)-user->Vp*user->tau_0*t)/user->lT;

				aF[j][i].U=(U1*(-0.5*(Ux*Px+Uy*Py)+0.5*(1-aY[j][i].p)*(grad))+U2*(1+(1-user->k)*aY[j][i].U)*(Pxt*Px+Pyt*Py+aYdot[j][i].p*(grad_2))+U2*(1-user->k)*aYdot[j][i].p*(Ux*Px+Uy*Py)+0.5*(1+(1-user->k)*aY[j][i].U)*aYdot[j][i].p-0.5*(1+user->k-(1-user->k)*aY[j][i].p)*aYdot[j][i].U-(1-aY[j][i].p)*0.5*user->V_nd/user->W*(0.5*(asol_old[j][i].vx+asol_old[j][i-1].vx)*Ux+0.5*(asol_old[j][i].vx+asol_old[j-1][i].vy)*Uy))*hx2;//

				aF[j][i].p=( e*e*(grad_2)+2.0*e*ep*(thetax*Px+thetay*Py)+(e*epp+ep*ep)*(thetay*Px-thetax*Py)+(aY[j][i].p-aY[j][i].p*aY[j][i].p*aY[j][i].p)-user->lambda*(aY[j][i].p*aY[j][i].p-1)*(aY[j][i].p*aY[j][i].p-1)*(aY[j][i].U+P1+0.01*arand[j][i]*aY[j][i].p)-(1)*(1.0+user->epsilon_tau*cos(4.0*theta))*aYdot[j][i].p )*hx2 ;


				if (asol_old[j][i].p<-0.9) {

				if (t<3.0){ 
//			ierr= PetscPrintf(PETSC_COMM_WORLD,"first if\t") ;CHKERRQ(ierr);
						aF[j][i].pp=aY[j][i].pp;aF[j][i].vx=aY[j][i].vx;aF[j][i].vy=aY[j][i].vy;}
					else {
//			ierr= PetscPrintf(PETSC_COMM_WORLD,"second if!\t") ;CHKERRQ(ierr);
					RR=i+2;LL=i-2;
					TT=j+2;BB=j-2;
					if (i==1) {LL=info2->mx-2;}//for moment x ll is info2->mx-2 but for poisson ll is info2->mx-1  
					if (i==info2->mx-2) {RR=0;} //for poisson rr is 0
					if (j==1) {BB=1;} //for moment y, bb is 1 but for poisson bb is 0 
					if (j==info2->my-2) {TT=info2->my-1;} //for poisson tt is info2->my-1

					Px_old=(asol_old[j][i].p-asol_old[j][i-1].p)/(user->hx);
					Py_old=(0.5*(asol_old[j+1][i-1].p+asol_old[j+1][i].p)-0.5*(asol_old[j-1][i-1].p+asol_old[j-1][i].p))/(2.0*user->hy);					
					pp_x=(aY[j][i].pp-aY[j][i-1].pp)/(user->hx);

					VVx_x=(aY[j][i+1].vx-aY[j][i-1].vx)/(2.0*user->hx);
					VVx_y=(aY[j+1][i].vx-aY[j-1][i].vx)/(2.0*user->hy);
					VVx_xx=((aY[j][i+1].vx+aY[j][i-1].vx-2.0*aY[j][i].vx)/(hx2));
					VVx_yy=((aY[j+1][i].vx+aY[j-1][i].vx-2.0*aY[j][i].vx)/hy2);
					grad_VVx=VVx_xx+VVx_yy;//(4.0*(aY[j][i+1].vx+aY[j][i-1].vx+aY[j+1][i].vx+aY[j-1][i].vx)+1.0*(aY[j+1][i+1].vx+aY[j+1][i-1].vx+aY[j-1][i+1].vx+aY[j-1][i-1].vx)-20.0*aY[j][i].vx)/(6.0*hx2);

					grad_2=(0.5*(asol_old[j][i+1].p+asol_old[j][i].p)+0.5*(asol_old[j][i-1].p+asol_old[j][LL].p)-2*0.5*(asol_old[j][i-1].p+asol_old[j][i].p))/hx2+(0.5*(asol_old[j+1][i-1].p+asol_old[j+1][i].p)+0.5*(asol_old[j-1][i-1].p+asol_old[j-1][i].p)-2*0.5*(asol_old[j][i].p+asol_old[j][i-1].p))/hy2;
					x1=0.5*aYdot[j][i].vx*(1-0.5*(asol_old[j][i].p+asol_old[j][i-1].p));
					x2=-0.5*aY[j][i].vx*0.5*(aYdot[j][i].p+aYdot[j][i-1].p);
					x3=-0.5*(1-0.5*(asol_old[j][i].p+asol_old[j][i-1].p))*pp_x;
					x4=user->nu_nd*0.5*(1-0.5*(asol_old[j][i].p+asol_old[j][i-1].p))*grad_VVx;
					x5=-user->nu_nd*(Px_old*VVx_x+Py_old*VVx_y);
					x6=-user->nu_nd*0.5*(aY[j][i].vx)*grad_2;
					x7=-0.25*user->h*user->nu_nd*aY[j][i].vx*(1-0.5*(asol_old[j][i].p+asol_old[j][i-1].p))*(1+0.5*(asol_old[j][i].p+asol_old[j][i-1].p))*(1+0.5*(asol_old[j][i].p+asol_old[j][i-1].p));
					x8=-0.5*(1-0.5*(asol_old[j][i].p+asol_old[j][i-1].p))*(aY[j][i].vx*VVx_x+0.25*(aY[j][i].vy+aY[j+1][i].vy+aY[j][i-1].vy+aY[j+1][i-1].vy)*VVx_y);
					x9=0.0;
					aF[j][i].vx=(x3+x4+x5+x6+x7+x8+x9-x1-x2)*user->hx;

					Px_old=(0.5*(asol_old[j][i+1].p+asol_old[j-1][i+1].p)-0.5*(asol_old[j][i-1].p+asol_old[j-1][i-1].p))/(2.0*user->hy);					
					Py_old=(asol_old[j][i].p-asol_old[j-1][i].p)/(user->hy);
					pp_y=(aY[j][i].pp-aY[j-1][i].pp)/(user->hy);

					VVy_x=(aY[j][i+1].vy-aY[j][i-1].vy)/(2.0*user->hx);
					VVy_y=(aY[j+1][i].vy-aY[j-1][i].vy)/(2.0*user->hy);
					VVy_xx=((aY[j][i+1].vy+aY[j][i-1].vy-2.0*aY[j][i].vy)/(hx2));
					VVy_yy=((aY[j+1][i].vy+aY[j-1][i].vy-2.0*aY[j][i].vy)/hy2);
					grad_VVy=VVy_xx+VVy_yy;//(4.0*(aY[j][i+1].vy+aY[j][i-1].vy+aY[j+1][i].vy+aY[j-1][i].vy)+1.0*(aY[j+1][i+1].vy+aY[j+1][i-1].vy+aY[j-1][i+1].vy+aY[j-1][i-1].vy)-20.0*aY[j][i].vy)/(6.0*hx2);
					grad_2=(0.5*(asol_old[j+1][i].p+asol_old[j][i].p)+0.5*(asol_old[j-1][i].p+asol_old[BB][i].p)-2*0.5*(asol_old[j][i-1].p+asol_old[j][i].p))/hx2+(0.5*(asol_old[j][i+1].p+asol_old[j-1][i+1].p)+0.5*(asol_old[j][i-1].p+asol_old[j-1][i-1].p)-2*0.5*(asol_old[j-1][i].p+asol_old[j][i].p))/hy2;
					y1=0.5*aYdot[j][i].vy*(1-0.5*(asol_old[j][i].p+asol_old[j-1][i].p));
					y2=-0.5*aY[j][i].vy*0.5*(aYdot[j][i].p+aYdot[j-1][i].p);
					y3=-0.5*(1-0.5*(asol_old[j][i].p+asol_old[j-1][i].p))*pp_y;
					y4=user->nu_nd*0.5*(1-0.5*(asol_old[j][i].p+asol_old[j-1][i].p))*grad_VVy;
					y5=-user->nu_nd*(Px_old*VVy_x+Py_old*VVy_y);
					y6=-user->nu_nd*0.5*(aY[j][i].vy)*grad_2;
					y7=-0.25*user->h*user->nu_nd*aY[j][i].vy*(1-0.5*(asol_old[j][i].p+asol_old[j-1][i].p))*(1+0.5*(asol_old[j][i].p+asol_old[j-1][i].p))*(1+0.5*(asol_old[j][i].p+asol_old[j-1][i].p));
					y8=-0.5*(1-0.5*(asol_old[j][i].p+asol_old[j-1][i].p))*(0.25*(aY[j][i].vx+aY[j][i+1].vx+aY[j-1][i].vx+aY[j-1][i+1].vx)*VVy_x+aY[j][i].vy*VVy_y);
					y9=-0.5*(1-0.5*(asol_old[j][i].p+asol_old[j-1][i].p))*user->Density_var*user->accelerate*user->cl0*((1+(1-user->k)*0.5*(asol_old[j][i].U+asol_old[j-1][i].U))*(1+user->k-(1-user->k)*0.5*(asol_old[j][i].p+asol_old[j-1][i].p))/(2.0*user->k)-1);
					aF[j][i].vy=(y3+y4+y5+y6+y7+y8+y9-y1-y2)*user->hx;//1- -((1+(1-user->k)*user->U_l)*2*user->cl0/(2*user->k))



//                                if (t<3.0){
//                                               aF[j][i].pp=aY[j][i].pp;}
//                                        else {
					if (i==1) {LL=info2->mx-1;}//for moment x ll is info2->mx-2 but for poisson ll is info2->mx-1
					if (j==1) {BB=0;} 
					Px_old=(asol_old[j][i+1].p-asol_old[j][i-1].p)/(2.0*user->hx);
					Py_old=(asol_old[j+1][i].p-asol_old[j-1][i].p)/(2.0*user->hy);
					Pxx_old=((asol_old[j][i+1].p+asol_old[j][i-1].p-2.0*asol_old[j][i].p)/(hx2));
					Pxy_old=((asol_old[j+1][i+1].p+asol_old[j-1][i-1].p-asol_old[j+1][i-1].p-asol_old[j-1][i+1].p)/(4.0*hxy));
					Pyy_old=((asol_old[j+1][i].p+asol_old[j-1][i].p-2.0*asol_old[j][i].p)/hy2);
					Pxxx_old=(asol_old[j][RR].p-2.0*asol_old[j][i+1].p+2.0*asol_old[j][i-1].p-asol_old[j][LL].p)/(2.0*hx3);
					Pyyy_old=(asol_old[TT][i].p-2.0*asol_old[j+1][i].p+2.0*asol_old[j-1][i].p-asol_old[BB][i].p)/(2.0*hy3);

					Pxxy_old= ((asol_old[j+1][i+1].p+asol_old[j+1][i-1].p-2.0*asol_old[j+1][i].p)-(asol_old[j-1][i+1].p+asol_old[j-1][i-1].p-2.0*asol_old[j-1][i].p))/(2.0*hx2y);
					Pxyy_old=((asol_old[j+1][i+1].p+asol_old[j-1][i+1].p-2.0*asol_old[j][i+1].p)-(asol_old[j+1][i-1].p+asol_old[j-1][i-1].p-2.0*asol_old[j][i-1].p))/(2.0*hxy2);
					
					Uy_old=(asol_old[j+1][i].U-asol_old[j-1][i].U)/(2.0*user->hy);
					
					pp_x=(aY[j][i+1].pp-aY[j][i-1].pp)/(2.0*user->hx);
					pp_y=(aY[j+1][i].pp-aY[j-1][i].pp)/(2.0*user->hy);
					pp_xx=((aY[j][i+1].pp+aY[j][i-1].pp-2.0*aY[j][i].pp)/(hx2));
					pp_yy=((aY[j+1][i].pp+aY[j-1][i].pp-2.0*aY[j][i].pp)/hy2);
					VVx_x=(aY[j][i+1].vx-aY[j][i].vx)/(user->hx);
					VVx_y=(0.5*(aY[j+1][i+1].vx+aY[j+1][i].vx)-0.5*(aY[j-1][i+1].vx+aY[j-1][i].vx))/(2.0*user->hy);
					VVx_xx=((0.5*(aY[j][i+1].vx+aY[j][RR].vx)+0.5*(aY[j][i-1].vx+aY[j][i].vx)-2*0.5*(aY[j][i+1].vx+aY[j][i].vx))/(hx2));
					//VVx_xx=((0.5*(aY[j][i].vx+aY[j][i].vx)+0.5*(aY[j][i].vx+aY[j][i].vx)-2*0.5*(aY[j][i].vx+aY[j][i].vx))/(hx2));
					VVx_yy=((0.5*(aY[j+1][i].vx+aY[j+1][i+1].vx)+0.5*(aY[j-1][i].vx+aY[j-1][i+1].vx)-2*0.5*(aY[j][i].vx+aY[j][i+1].vx))/(hy2));
					grad_VVx=VVx_xx+VVx_yy;//(4.0*(aY[j][i+1].vx+aY[j][i-1].vx+aY[j+1][i].vx+aY[j-1][i].vx)+1.0*(aY[j+1][i+1].vx+aY[j+1][i-1].vx+aY[j-1][i+1].vx+aY[j-1][i-1].vx)-20.0*aY[j][i].vx)/(6.0*hx2);
					VVy_x=(0.5*(aY[j][i+1].vy+aY[j+1][i+1].vy)-0.5*(aY[j][i-1].vy+aY[j+1][i-1].vy))/(2.0*user->hx);//0.5*(aY[j][i].vy+aY[j][i].vy)
					VVy_y=(aY[j+1][i].vy-aY[j][i].vy)/(user->hy);
					VVy_xx=((0.5*(aY[j][i+1].vy+aY[j+1][i+1].vy)+0.5*(aY[j][i-1].vy+aY[j+1][i-1].vy)-2.0*0.5*(aY[j+1][i].vy+aY[j][i].vy))/(hx2));
					VVy_yy=((0.5*(aY[j+1][i].vy+aY[TT][i].vy)+0.5*(aY[j][i].vy+aY[j-1][i].vy)-2.0*0.5*(aY[j+1][i].vy+aY[j][i].vy))/hy2);
					grad_VVy=VVy_xx+VVy_yy;//(4.0*(aY[j][i+1].vy+aY[j][i-1].vy+aY[j+1][i].vy+aY[j-1][i].vy)+1.0*(aY[j+1][i+1].vy+aY[j+1][i-1].vy+aY[j-1][i+1].vy+aY[j-1][i-1].vy)-20.0*aY[j][i].vy)/(6.0*hx2);

//0.5*(aY[j][i].vx+aY[j][i+1].vx)      
/*					if (t>0.000001){ 
						a1=-0.5*(0.5*(aYdot[j][i].vx+aYdot[j][i+1].vx)*Px_old+0.5*(aYdot[j][i].vy+aYdot[j+1][i].vy)*Py_old);
						a9=-user->Density_var*user->accelerate*user->cl0*(  -0.5*Py_old*((1+(1-user->k)*asol_old[j][i].U)*(1+user->k-(1-user->k)*asol_old[j][i].p)/2.0/user->k -1)+0.5*(1-asol_old[j][i].p)*((1-user->k)/(2.0*user->k))*Uy_old*(1+user->k-(1-user->k)*asol_old[j][i].p)-0.5*(1-asol_old[j][i].p)*((1-user->k)/(2.0*user->k))*Py_old*(1+(1-user->k)*asol_old[j][i].U));
						}
					else {
						a1=0.0;a9=0.0;
						}*/


					a1=-0.5*(0.5*(aYdot[j][i].vx+aYdot[j][i+1].vx)*Px_old+0.5*(aYdot[j][i].vy+aYdot[j+1][i].vy)*Py_old);
					a9=-user->Density_var*user->accelerate*user->cl0*(  -0.5*Py_old*((1+(1-user->k)*asol_old[j][i].U)*(1+user->k-(1-user->k)*asol_old[j][i].p)/2.0/user->k -1)+0.5*(1-asol_old[j][i].p)*((1-user->k)/(2.0*user->k))*Uy_old*(1+user->k-(1-user->k)*asol_old[j][i].p)-0.5*(1-asol_old[j][i].p)*((1-user->k)/(2.0*user->k))*Py_old*(1+(1-user->k)*asol_old[j][i].U));
					a2=-0.5*(0.5*(aY[j][i].vx+aY[j][i+1].vx)*Pxt+0.5*(aY[j][i].vy+aY[j+1][i].vy)*Pyt);
					a3=0.5*(1-asol_old[j][i].p)*(2*VVx_y*VVy_x-2.0*VVx_x*VVy_y);////////
					a4=-0.5*Px_old*(0.5*(aY[j][i].vx+aY[j][i+1].vx)*VVx_x+0.5*(aY[j][i].vy+aY[j+1][i].vy)*VVx_y)-0.5*Py_old*(0.5*(aY[j][i].vx+aY[j][i+1].vx)*VVy_x+0.5*(aY[j][i].vy+aY[j+1][i].vy)*VVy_y);
					a5=-0.5*(1-asol_old[j][i].p)*(pp_xx+pp_yy);
					a6=0.5*(pp_x*Px_old+pp_y*Py_old);
					a7=-(0.5*Px_old*grad_VVx+0.5*Py_old*grad_VVy+Pxy_old*(VVx_y+VVy_x)+Pxx_old*VVx_x+Pyy_old*VVy_y+0.5*0.5*(aY[j][i].vy+aY[j+1][i].vy)*(Pxxy_old+Pyyy_old)+0.5*0.5*(aY[j][i].vx+aY[j][i+1].vx)*(Pxyy_old+Pxxx_old))*user->nu_nd;
					a8=-0.25*user->h*user->nu_nd*(1+asol_old[j][i].p)*(1-3.0*asol_old[j][i].p)*(0.5*(aY[j][i].vx+aY[j][i+1].vx)*Px_old+0.5*(aY[j][i].vy+aY[j+1][i].vy)*Py_old);
//replace 0.00 by  
					aF[j][i].pp=(a5+a6+a7+a8+a9-a1-a2-a3-a4)*hx2;}// replace 0.0 with a9
					}
					else {

					aF[j][i].vx=aY[j][i].vx;
					aF[j][i].vy=aY[j][i].vy;
					aF[j][i].pp=aY[j][i].pp;
						}

						}
//			ierr= PetscPrintf(PETSC_COMM_WORLD,"end (%d,%d)!\n",i,j) ;CHKERRQ(ierr);
					}
				}
//			}


	ierr= DMDAVecRestoreArray(user->da,user->rand,&arand) ;CHKERRQ(ierr);
	ierr= DMDAVecRestoreArray(user->da,user->sol_local,&asol_old) ;CHKERRQ(ierr);

//	ierr= PetscPrintf(PETSC_COMM_WORLD,"end time= %f!\n",t) ;CHKERRQ(ierr);
	return 0;
	}

int main (int argc, char **argv){


//variable definition
PetscErrorCode ierr;
DMDALocalInfo info2;

//initialize petsc
ierr= PetscInitialize(&argc,&argv, (char*) 0,help);CHKERRQ(ierr);

//variable initialization
struct VAR_STRUCT user;
var_input(&user);

//mpi_size


ierr= PetscSNPrintf(user.fname,sizeof(user.fname),"result.dat");CHKERRQ(ierr);
ierr= PetscFOpen(PETSC_COMM_WORLD,user.fname,"a",&user.fp);CHKERRQ(ierr);
ierr= PetscFPrintf(PETSC_COMM_WORLD,user.fp," Time[DL]\tn\tTime step[DL]\tsolution-step[DL]\tTS-step\n");



//DMDA Create
ierr= DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,
                          user.Nx,user.Ny,PETSC_DECIDE,PETSC_DECIDE,5,4,NULL,NULL,&user.da) ;CHKERRQ(ierr);    //
ierr= DMSetUp(user.da) ;CHKERRQ(ierr);
ierr= DMDASetFieldName(user.da,0,"U") ;CHKERRQ(ierr); //set the name of the first field component temperature
ierr= DMDASetFieldName(user.da,1,"p") ;CHKERRQ(ierr); //set the name of the second field component order parameter
ierr= DMDASetFieldName(user.da,2,"vy") ;CHKERRQ(ierr); //set the name of the second field component order parameter
ierr= DMDASetFieldName(user.da,3,"vx") ;CHKERRQ(ierr); //set the name of the second field component order parameter
ierr= DMDASetFieldName(user.da,4,"pp") ;CHKERRQ(ierr); //set the name of the second field component order parameter

ierr= DMSetApplicationContext(user.da,&user) ;CHKERRQ(ierr); //set  a user context to a DM object
ierr= DMDAGetLocalInfo(user.da,&info2) ;CHKERRQ(ierr);

//vec create
ierr= DMCreateGlobalVector(user.da,&user.sol) ;CHKERRQ(ierr);
ierr= DMCreateGlobalVector(user.da,&user.sol_old) ;CHKERRQ(ierr);
ierr= DMCreateLocalVector(user.da,&user.sol_local) ;CHKERRQ(ierr);
ierr= DMCreateGlobalVector(user.da,&user.rand) ;CHKERRQ(ierr);

//random number 
ierr=PetscRandomCreate(PETSC_COMM_WORLD,&user.r);CHKERRQ(ierr);
ierr=PetscRandomSetType(user.r,PETSCRAND48);CHKERRQ(ierr);
ierr=PetscRandomSetInterval(user.r,-0.5,0.5);CHKERRQ(ierr);


//TS 
ierr= TSCreate(PETSC_COMM_WORLD, &user.ts) ; CHKERRQ(ierr);
ierr= TSSetProblemType(user.ts,TS_NONLINEAR) ; CHKERRQ(ierr);
ierr= TSSetDM(user.ts,user.da); CHKERRQ(ierr);
ierr= DMDATSSetIFunctionLocal(user.da,INSERT_VALUES,(DMDATSIFunctionLocal)FormFunction,&user) ; CHKERRQ(ierr);
ierr= TSSetExactFinalTime(user.ts,TS_EXACTFINALTIME_MATCHSTEP) ; CHKERRQ(ierr);
ierr= TSSetTimeStep(user.ts,user.ht) ; CHKERRQ(ierr);
//ierr= TSSetMaxSteps(user.ts,200000) ; CHKERRQ(ierr);
ierr= TSSetDuration(user.ts,user.Nt,user.t_tot) ; CHKERRQ(ierr);
ierr= TSSetFromOptions(user.ts) ; CHKERRQ(ierr);
//initialize vectors
ierr= InitialState(user.sol,&info2,&user) ;CHKERRQ(ierr); //initialize previously
ierr= VecCopy(user.sol,user.sol_old); CHKERRQ(ierr);

//solve

ierr = TSSolve(user.ts,user.sol); CHKERRQ(ierr);

//destroy
ierr= VecDestroy(&user.sol);CHKERRQ(ierr);
ierr= VecDestroy(&user.sol_old);CHKERRQ(ierr);
ierr= VecDestroy(&user.rand);CHKERRQ(ierr);
ierr= PetscFinalize();CHKERRQ(ierr);
return 0;
}