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


		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);CHKMEMQ;

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

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

	hx2=user->hx*user->hx;CHKMEMQ;
	hy2=user->hy*user->hy;CHKMEMQ;
	hxy=user->hx*user->hy;CHKMEMQ;
	hx3=user->hx*user->hx*user->hx;CHKMEMQ;
	hy3=user->hy*user->hy*user->hy;CHKMEMQ;
	hx2y=user->hx*user->hx*user->hy;CHKMEMQ;
	hxy2=user->hx*user->hy*user->hy;CHKMEMQ;
	ierr= DMDAVecGetArray(user->da,user->rand,&arand) ;CHKERRQ(ierr);CHKMEMQ;
	ierr= DMDAVecGetArray(user->da,user->sol_local,&asol_old) ;CHKERRQ(ierr);CHKMEMQ;
	for (j=info2->ys;j<info2->ys+info2->ym;j++){
		for (i=info2->xs;i<info2->xs+info2->xm;i++){
			 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;CHKMEMQ;}
			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;CHKMEMQ;}
			else {
				aF[j][i].U=aY[j][i].U;
				aF[j][i].p=aY[j][i].p ;
				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= DMDAVecRestoreArray(user->da,user->rand,&arand) ;CHKERRQ(ierr);
	ierr= DMDAVecRestoreArray(user->da,user->sol_local,&asol_old) ;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,3,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;
}