static char help[] = "Solving swift hohenberg equation in 1D \n\
u_t=(r-1)u-2u_xx+u_xxxx in 2D.\n";

#include <petscts.h>
#include <petscdm.h>
#include <petscdmda.h>

extern PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void*);
extern PetscErrorCode InitialValues(DM, Vec, void*);
extern PetscErrorCode RHSJacobian(TS, PetscReal, Vec, Mat, Mat, void*);
extern PetscErrorCode Monitor(TS,PetscInt,PetscReal,Vec,void*);

int main(int argc, char **argv)
{
  Vec x;
  Mat J;
  TS ts;
  PetscReal a,b,c,h,ftime, dt;
  DM da;
  PetscErrorCode ierr;
  PetscInt  steps, maxsteps=100;

  PetscInitialize(&argc, &argv, (char*)0,help);

  /*==================================================================	
    Create 1D distributed array for managing parallel grid and vectors	
    ==================================================================*/

  ierr = DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 15, 1, 1, NULL, &da); CHKERRQ(ierr);

  /*==================================================================	
    Let DMDA create vector						
    ==================================================================*/

  ierr = DMCreateGlobalVector(da, &x);CHKERRQ(ierr);

  /*==================================================================	
    Create timestepper context to use Euler Method		
    ==================================================================*/

  ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
  ierr = TSSetDM(ts, da);CHKERRQ(ierr);
  ierr = TSSetType(ts, TSEULER);CHKERRQ(ierr);

  ierr = TSMonitorSet(ts, Monitor, NULL, NULL);

  ierr = TSSetRHSFunction(ts, x, RHSFunction, NULL);CHKERRQ(ierr);

  /*=================================================================
    Create Jacobian for Euler Method
    =================================================================*/
  ierr = DMSetMatType(da, MATAIJ);CHKERRQ(ierr);
  ierr = DMCreateMatrix(da,&J);CHKERRQ(ierr);
  ierr = TSSetRHSJacobian(ts,J,J,RHSJacobian,NULL); CHKERRQ(ierr);

  ftime = 1.0;
  ierr = TSSetDuration(ts, maxsteps, ftime);CHKERRQ(ierr);
  ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr);

  /*==================================================================	
    Set initial conditions (see user function InitialValues();)		
    ==================================================================*/

  ierr = InitialValues(da, x, NULL);CHKERRQ(ierr);
  // ierr = VecView(x, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  dt = 0.0001;
  ierr = TSSetInitialTimeStep(ts, 0.0, dt); CHKERRQ(ierr);

  /*==================================================================	
    Set runtime options for TS
    ==================================================================*/
  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);

  //ierr = MatView(J, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

  /*==================================================================	
    Solve PDE at each time step
    ==================================================================*/
  ierr = PetscPrintf(PETSC_COMM_WORLD, "time stepping\n\n");CHKERRQ(ierr);
  ierr = TSSolve(ts, x);CHKERRQ(ierr);
  ierr = TSGetSolveTime(ts, &ftime);CHKERRQ(ierr);
  ierr = TSGetTimeStepNumber(ts, &steps);CHKERRQ(ierr);


   ierr = MatView(J, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  /*==================================================================	
    Clean up created objects
    ==================================================================*/

  ierr = VecDestroy(&x);CHKERRQ(ierr);
  ierr = MatDestroy(&J);CHKERRQ(ierr);
  ierr = TSDestroy(&ts);CHKERRQ(ierr);
  ierr = DMDestroy(&da);CHKERRQ(ierr);

  PetscFinalize();
  PetscFunctionReturn(0);
}


  /*==================================================================	
    Function to evaluate RHS of each time step
    ==================================================================*/
#undef __FUNCT__	       
#define __FUNCT__ "RHSFunction"

PetscErrorCode RHSFunction(TS ts, PetscReal ftime, Vec X, Vec F, void *unused)
{
  DM da;
  PetscErrorCode ierr;
  PetscInt i, Mx, xs, xm;
  PetscReal h;
  PetscScalar u, uxx, uxxxx,  *uarray, *f;
  Vec localX;
  
  /* Always begin user defined functions with this call */
  PetscFunctionBeginUser;

  /* Access DM for TS  */
  ierr = TSGetDM(ts, &da); CHKERRQ(ierr);
  
  /* Get a local part of the vecor x  */
  ierr = DMGetLocalVector(da, &localX);CHKERRQ(ierr);

  /* Access info about da */
  ierr = DMDAGetInfo(da, PETSC_IGNORE, &Mx, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE); CHKERRQ(ierr); 

  /* Scatter ghost points */
  ierr = DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX);CHKERRQ(ierr);
  ierr = DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX);CHKERRQ(ierr);

  /* Give local DM access to local vector data */
  ierr = DMDAVecGetArrayRead(da, localX, &uarray);CHKERRQ(ierr);
  ierr = DMDAVecGetArray(da, F, &f);CHKERRQ(ierr);

  /* Get where the local grid boundaries are */
  ierr = DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL);CHKERRQ(ierr);

  /* Compute point separation in discretization */
  h = 1/(Mx-1);
  
  //  ierr = PetscPrintf(PETSC_COMM_WORLD, "Calculating RHS Functino\n\n");CHKERRQ(ierr);
  /* Compute the function over the local grid points */
  for (i=xs; i<xs+xm; i++) {
    if(i==0 || i==Mx-1) {
      f[i]=uarray[i];
      continue;
    }
    u = uarray[i];
    uxx = (uarray[i+1]-2.0*uarray[i]+uarray[i-1])/(h*h);
    uxxxx = (uarray[i+2]-4.0*uarray[i+1]+6.0*uarray[i]-4.0*uarray[i-1]+uarray[i-2])/(h*h*h*h);
    f[i] =u-2.0*uxx+uxxxx;
  }
  
  /* Restore the vectors */
  ierr = DMDAVecRestoreArrayRead(da,localX, &uarray);CHKERRQ(ierr);
  ierr = DMDAVecRestoreArray(da, F, &f);CHKERRQ(ierr);
  ierr = DMRestoreLocalVector(da, &localX);CHKERRQ(ierr);

  /* end user defined function */
  PetscFunctionReturn(0);
}

  /*==================================================================	
    User defined function to define initial values
    ==================================================================*/
#undef __FUNCT__	       
#define __FUNCT__ "InitialValues"

PetscErrorCode InitialValues(DM da, Vec X, void *unused)
{
  PetscErrorCode ierr;
  PetscInt i, xs, xm, Mx;
  PetscScalar *x;
  PetscReal h;
  /* Always start user defined functions with this call */
  PetscFunctionBeginUser;

  /* Get DMDA info for global and local vector lengths */
  ierr = DMDAGetInfo(da, PETSC_IGNORE, &Mx, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE);

  /* Get vector data for local operations  */
  ierr = DMDAVecGetArray(da, X, &x); CHKERRQ(ierr);

  /* Get corner data to know length of local vector */
  ierr = DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL); CHKERRQ(ierr);
  
  /* Set initial values */
  h = 1./(Mx-1);
  for (i=xs; i<xs+xm; i++) {
    x[i] =100.*sin(2*3.1415*i*h);
  }

  /* Restore array to update values */
  ierr = DMDAVecRestoreArray(da, X, &x);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

  /*==================================================================	
    User provided routine for RHS's Jacobian
    ==================================================================*/
#undef __FUNCT__	       
#define __FUNCT__ "RHSJacobian"

PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec X, Mat J, Mat Jpre, void *unused)
{
  PetscErrorCode ierr;
  DM da;
  DMDALocalInfo info;
  PetscInt i;
  PetscReal h;
  
  
  /* Always start user defined functions with this call */
  PetscFunctionBeginUser;

  /* Get TS info and DMDA info */
  ierr = TSGetDM(ts, &da);CHKERRQ(ierr);
  ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);

  /* Set local grid spacing */
  h = 1.0/(info.mx-1);
  PetscPrintf(PETSC_COMM_WORLD,"Definining Jacobian");
  /*for loop for assigning matrix values */
  for (i=info.xs; i<info.xs+info.xm;i++) {
    PetscInt nc=0;
    MatStencil row,col[5];
    PetscScalar val[5];
    row.i=i;
    if (i==0 || i==info.mx-1) {
      col[nc].i=i; val[nc++]=1.0;
    } else {
      col[nc].i=i-2;  val[nc++]=1./(h*h*h*h);
      col[nc].i=i-1;  val[nc++]=-2./(h*h)-4./(h*h*h*h);
      col[nc].i=i;    val[nc++]=1+4./(h*h)+6./(h*h*h*h);
      col[nc].i=i+1;  val[nc++]=-2./(h*h)-4./(h*h*h*h);
      col[nc].i=i+2;  val[nc++]=1./(h*h*h*h);
    }
    ierr = MatSetValuesStencil(Jpre, 1, &row, nc, col, val, INSERT_VALUES);CHKERRQ(ierr);
  }

  ierr = MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  if (J!=Jpre) {
    ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    ierr = MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);  
  }
  PetscFunctionReturn(0);
}

#undef __FUNCT__	       
#define __FUNCT__ "Monitor"

PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal time, Vec X, void* unused){
  PetscFunctionBeginUser;
  VecView(X, PETSC_VIEWER_DRAW_WORLD);
  PetscFunctionReturn(0);
}