#include <petscdm.h>
#include <petscdmlabel.h>
#include <petscds.h>
#include <petscdmplex.h>
#include <petscksp.h>
#include <petscsnes.h>
#include <petscts.h>


typedef enum {RUN_STANDALONE, RUN_COUPLED} RunType;
typedef enum {NONE, LABELS, BOUNDARIES, SOLUTION} Visualization;


typedef struct  /* The time context */
{
  PetscReal total;  /* total time in seconds */
  PetscReal dt;                      /* the current time step size */
  PetscInt iter;                     /* number of iterations so far */
  /* PetscBool halfstep;  /\* Is the timestep a halfstep? *\/ */
} TSctx;


typedef struct {
  PetscInt debug;
  RunType runType;
  /* definitions for mesh and problem setup */
  PetscBool simplex;
  PetscBool interpolate;
  /* Additional options and parameters */
  PetscBool verbose;
  TSctx time;
  PetscViewer viewer;
}AppCtx;


/* defining the kronecker delta */
const PetscInt delta3D[3*3] = {1,0,0,0,1,0,0,0,1};


/* defining the material parameters */
const PetscReal mu =76.923076923, lbda=115.384615385, rho_steel = 7850E-9;



void f0_u_transient(PetscInt dim, PetscInt Nf, PetscInt NfAux,
                    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[],
                    const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[],
                    const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[],
                    const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscScalar f0[])
{
  const PetscInt Ncomp = dim;
  PetscInt comp;
  for(comp=0;comp<Ncomp;comp++) f0[comp]= 1.0*rho_steel*u_t[uOff[1]+comp];
}



void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
          const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[],
          const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[],
          const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[],
          const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscScalar f1[])
{
  const PetscInt Ncomp = dim;
  
  /* f1 is the cauchy stress tensor*/
  
  PetscInt comp, d;
  for(comp=0;comp<Ncomp;comp++){
    for(d=0;d<dim;d++){
      f1[comp*dim+d]=mu*(u_x[comp*dim+d]+u_x[d*dim+comp]);
    }
    for(d=0;d<dim;d++){
      f1[comp*dim+comp]+=lbda*u_x[d*dim+d];
    }
  }
}

void g3_uu_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux,
	      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], 
	      const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[],
	      const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[],
	      const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[],
	      PetscScalar g3[]){
  
  const PetscInt Ncomp = dim;
   PetscInt i,j,k,l,m,n;
    PetscScalar G[dim*dim*Ncomp*Ncomp];  /* Cauchy-Green  strain tensor */
  PetscScalar D[dim*dim*Ncomp*Ncomp];  /* constitutive tensor for isotropic material */
    

  /* partial derivative of the cauchy strain ensor with respect to the displacement gradient G_mjnl = (partial epsilon_mn)/(partial u^j_l) */
  for (m=0;m<dim;m++){  /* primary index of the cauchy green strain tensor */
    for (j=0;j<dim;j++){  /* component of trial function */
      for (n=0;n<Ncomp;n++){  /* secondary index of the cauchy green strain tensor */
        for (l=0;l<Ncomp;l++){  /* derivative index for trial function */
          G[((m*dim+j)*dim+n)*dim+l]=0.5*(delta3D[n*dim+l]*delta3D[j*dim+m]+delta3D[m*dim+l]*delta3D[n*dim+j]);
        }
      }
    }
  }

  /* constructing the constitutive tensor for isotropic material */
  for (i=0;i<dim;i++){  /* component of test function */
    for (m=0;m<dim;m++){  /* primary index of the cauchy strain tensor */
      for (k=0;k<dim;k++){  /* derivative index of the test function */
        for (n=0;n<dim;n++){  /* secondary index of the cauchy strain
                               * tensor*/
          D[((i*dim+m)*dim+k)*dim+n] = lbda*delta3D[i*dim+k]*delta3D[m*dim+n]+mu*(delta3D[i*dim+m]*delta3D[k*dim+n]+delta3D[i*dim+n]*delta3D[k*dim+m]);
        }
      }
    }
  }
  /* constructing the integrand for gradient of testfunction and gradient of trialfunction */
  for (i=0;i<dim;i++){  /* component of test function */
    for (j=0;j<dim;j++){  /* component of the trial function */
      for (k=0;k<dim;k++){  /* derivative index of the test funciton */
        for (l=0;l<dim;l++){  /* derivative index of the trial function */
          for (m=0;m<dim;m++){  /* primary index of the cauchy strain
                                 * tensor */
            for (n=0;n<dim;n++){  /* secondary index of the cauchy
                                   * strain tensor */
              g3[((i*dim+j)*dim+k)*dim+l] +=  D[((i*dim+m)*dim+k)*dim+n]*(G[((m*dim+j)*dim+n)*dim+l]);
            }
          }
        }
      }
    }
  }
}


void f0_u_bd(PetscInt dim, PetscInt Nf, PetscInt NfAux,
                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[],
                const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[],
                const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[],
                const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[],
                PetscScalar f0[])
{
  const PetscInt Ncomp=dim;

  PetscInt comp;
  PetscReal pressure = 0.001;  /* in GPa */

  for (comp=0; comp<Ncomp; ++comp){
    f0[comp] = pressure*n[comp];
  }
  
}

    
void f1_u_bd(PetscInt dim, PetscInt Nf, PetscInt NfAux,
        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[],
        const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[],
        const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[],
        const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[],
        PetscScalar f1[])
    {
    const  PetscInt Ncomp = dim;
    PetscInt comp,d;
    
    for (comp=0;comp<Ncomp;++comp){
      for (d = 0; d<dim;++d){
        f1[comp*dim + d] = 0.0;
      }
    }
  }


void f0_vel(PetscInt dim, PetscInt Nf, PetscInt NfAux,
               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[],
               const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[],
               const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[],
               const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscScalar f0[])
{
  const PetscInt Ncomp = dim;
  PetscInt comp, d;
  for (comp=0;comp<Ncomp;comp++){
    f0[comp] = u[uOff[1]+comp]-u_t[comp];
  }
}


void f1_vel(PetscInt dim, PetscInt Nf, PetscInt NfAux,
               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[],
               const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[],
               const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[],
               const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscScalar f1[])
{
  const PetscInt Ncomp = dim;
  PetscInt comp,d;

  for (comp=0;comp<Ncomp;comp++){
    for (d=0;d<dim;d++){
      f1[comp*Ncomp+d] = 0.0;
    }
  }
}

void g0_velvel(PetscInt dim, PetscInt Nf, PetscInt NfAux,
                  const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], 
                  const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[],
                  const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[],
                  const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[],
                  PetscScalar g0[]){
  PetscInt Ncomp = dim;
  PetscInt i;

  for (i=0;i<Ncomp;i++){
    g0[i*Ncomp+i]= 1.0;
  }
  
}

void g0_uvel(PetscInt dim, PetscInt Nf, PetscInt NfAux,
                  const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], 
                  const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[],
                  const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[],
                  const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[],
                  PetscScalar g0[]){
  PetscInt Ncomp = dim;
  PetscInt i;

  for (i=0;i<Ncomp;i++){
    g0[i*Ncomp+i]= 1.0*rho_steel;
  }
  
}

void g2_velu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
                  const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], 
                  const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[],
                  const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[],
                  const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[],
                  PetscScalar g2[]){

  const PetscInt NComp=dim;
  PetscInt i,j,k;
  /* This is for the material derivative */
  for (i=0;i<dim;i++){
    for (j=0; j < dim; ++j) {
      for (k=0; k<dim; k++){
        g2[(i*dim+j)*dim+k] = -u_t[i];
      }
    }
  }

}


void g0_velu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
                  const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], 
                  const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[],
                  const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[],
                  const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[],
                  PetscScalar g0[]){

  PetscInt Ncomp = dim;
  PetscInt i,j;

  for (i=0;i<Ncomp;i++){
    g0[i*Ncomp+i]=-1.0;
  }
}

void f0_vel_bd(PetscInt dim, PetscInt Nf, PetscInt NfAux,
        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[],
        const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[],
        const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[],
        const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[],
        PetscScalar f0[])
{
  const  PetscInt Ncomp = dim;
  PetscInt comp,d;
  
  for (comp=0;comp<Ncomp;++comp){
    f0[comp] = 0.0;
  }
}

void f1_vel_bd(PetscInt dim, PetscInt Nf, PetscInt NfAux,
        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[],
        const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[],
        const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[],
        const PetscScalar a_x[], PetscReal t, const PetscReal x[], const PetscReal n[],
        PetscScalar f1[])
{
  const  PetscInt Ncomp = dim;
  PetscInt comp,d;
  
  for (comp=0;comp<Ncomp;++comp){
    for (d = 0; d<dim;++d){
      f1[comp*dim + d] = 0.0;
    }
  }
}



static PetscErrorCode zero_scalar(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
{
  u[0]=0.0;
  return 0;
}

static PetscErrorCode zero_vector(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
{
  const PetscInt Ncomp = dim;
  PetscInt comp;
  for (comp=0;comp<Ncomp;comp++) u[comp]=0.0;
  return 0;
}



static PetscErrorCode SetupProblem(DM dm, AppCtx *user)
{
  PetscDS prob; 		  
  PetscErrorCode ierr;
  PetscInt NumFields;


  ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);


  
     /* Setting up the linear elasticity problem */
  
  ierr = PetscDSSetResidual(prob,0,f0_u_transient,f1_u);CHKERRQ(ierr);
  ierr = PetscDSSetJacobian(prob,0,0,NULL,NULL,NULL,g3_uu_3d);CHKERRQ(ierr);

     


    /* setting up the velocity field for the transient analysis */
  
  ierr = PetscDSSetResidual(prob,1,f0_vel,f1_vel);CHKERRQ(ierr);
  ierr = PetscDSSetJacobian(prob,1,1,g0_velvel,NULL,NULL,NULL);CHKERRQ(ierr);
  ierr = PetscDSSetDynamicJacobian(prob,1,0,g0_velu,NULL,NULL,NULL);CHKERRQ(ierr);
  ierr = PetscDSSetDynamicJacobian(prob,0,1,g0_uvel,NULL,NULL,NULL);CHKERRQ(ierr);
  ierr = PetscDSSetBdResidual(prob,1,f0_vel_bd,f1_vel_bd);CHKERRQ(ierr);
  
  
  /* Setting the Neumann Boudnary Condition */
  
  ierr = PetscDSSetBdResidual(prob,0,f0_u_bd,f1_u_bd);CHKERRQ(ierr);

  
  ierr = PetscDSGetNumFields(prob,&NumFields);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"Problem set up with %i fields\n",NumFields);CHKERRQ(ierr);
      
  return(0);
}


static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user){

  DM cdm = dm;
  const PetscInt dim = 3;
  PetscFE fe, fe_bd, fe_vel, fe_vel_bd;
  PetscDS prob;
  PetscErrorCode ierr;
  PetscBool simplex = user->simplex;
  PetscQuadrature q;
  PetscInt order;
  PetscInt BSpaceOrder;
  PetscSpace space;

  /* Creating the FE */
  ierr = PetscFECreateDefault(dm, dim, dim, simplex,"def_",PETSC_DEFAULT,&fe);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) fe, "deformation");CHKERRQ(ierr);
  
  ierr = PetscFEGetQuadrature(fe,&q);CHKERRQ(ierr);
  ierr = PetscQuadratureGetOrder(q,&order);CHKERRQ(ierr);
  /* Creating the FE field for the velocity */
  ierr = PetscFECreateDefault(dm,dim,dim,simplex,"vel_",PETSC_DEFAULT,&fe_vel);CHKERRQ(ierr);
  ierr = PetscFESetQuadrature(fe_vel,q);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) fe_vel, "velocity");CHKERRQ(ierr);
    /* Discretization and boundary conditons: */
  while (cdm)  {
    ierr = DMGetDS(cdm, &prob);CHKERRQ(ierr);
    ierr = PetscDSSetDiscretization(prob, 0, (PetscObject) fe); CHKERRQ(ierr);
    ierr = PetscDSSetDiscretization(prob,1, (PetscObject) fe_vel);CHKERRQ(ierr);
    
    ierr = SetupProblem(cdm, user); CHKERRQ(ierr);
    /* Define the boundaries */
    
    const PetscInt Ncomp = dim;
    PetscInt components[dim];
    const PetscInt Nfid = 2;
    PetscInt d;
    PetscInt fid[Nfid];  /* fixed faces [numer of fixed faces] */
    const PetscInt Npid = 1;  /* number of pressure loaded faces */
    PetscInt pid[Npid];       /* the ids of the pressure loaded faces */
    const PetscInt Nsxid = 1; /* number of face ids for the symmetry in x
                               * direction bc */
    PetscInt sxid[Nsxid];
    const PetscInt Nsyid = 1;
    PetscInt syid[Nsyid];

    PetscInt restrictX[] = {0};  /* array of constricted components for the
                                  * dmaddboundary routine */
    PetscInt restrictY[] = {1};
    PetscInt restrictZ[] = {2}; 
    

    PetscInt restrictall[3] = {0, 1, 2};  /* restricting all movement */
    
    for (d=0;d<dim;d++){
      components[d]=d;
    }
    pid[0]= 3;  /* The pressure loaded faces */
    sxid[0] = 5;  /* symmetry in X-direction */
    fid[0] = 1;   /* fixed faces */
    fid[1] = 2;   /* fixed faces */
    syid[0] = 4;  /* symmetry in Y-direction */

    ierr = DMAddBoundary(cdm, PETSC_TRUE, "symmx", "Face Sets",0,1,restrictX, (void (*)()) zero_scalar, Nsxid, sxid, user);CHKERRQ(ierr);
    ierr = DMAddBoundary(cdm, PETSC_TRUE, "symmy", "Face Sets",0,1,restrictY, (void (*)()) zero_scalar, Nsyid, syid, user);CHKERRQ(ierr);
    ierr = DMAddBoundary(cdm, PETSC_TRUE, "bottom and top", "Face Sets",0,Ncomp,restrictall, (void (*)()) zero_vector, Nfid, fid, user);CHKERRQ(ierr);

    ierr = DMAddBoundary(cdm, PETSC_FALSE, "load", "Face Sets",0, Ncomp, components, NULL, Npid, pid, user);CHKERRQ(ierr);
    
    ierr = DMGetCoarseDM(cdm, &cdm); CHKERRQ(ierr);
  }

  ierr = PetscFEGetBasisSpace(fe,&space);CHKERRQ(ierr);
  ierr = PetscSpaceGetOrder(space,&BSpaceOrder);CHKERRQ(ierr);
  PetscPrintf(PETSC_COMM_WORLD,"The order of the basis space is %i \n", BSpaceOrder);

  ierr = PetscFEDestroy(&fe); CHKERRQ(ierr);
  
  ierr = PetscFEDestroy(&fe_vel); CHKERRQ(ierr);
  
  return(0);
}

static PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal crtime, Vec u, void *ctx)
{
  AppCtx *user=(AppCtx*) ctx;
  PetscErrorCode ierr;
  PetscViewer viewer;
  char filename[50];
 
  sprintf(filename,"solution%i.vtu",step);
  ierr = PetscViewerVTKOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
  if (user->verbose==PETSC_TRUE) PetscPrintf(PETSC_COMM_WORLD,"Writing the solution to the vtk file. \n");
  ierr = VecView(u,viewer);CHKERRQ(ierr);
  ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
  if (user->verbose==PETSC_TRUE)PetscPrintf(PETSC_COMM_WORLD,"Done creating the VTK file. \n");
    
  return(0);
}


static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
{
  const char *runTypes[2]={"standalone","coupled"};
  const char *visutype[4]={"none","labels","boundaries","solution"};
  PetscInt vis;
  PetscInt run;
  PetscErrorCode ierr;
  options->debug = 0;
  options->runType = RUN_STANDALONE;
  options->interpolate = PETSC_TRUE;






  options->verbose = PETSC_FALSE;




  options->time.total = 0.001;
  options->time.dt = 0.000005;

   ierr = PetscOptionsBegin(comm,"", "Linear elasticity problem options", "DMPLEX"); CHKERRQ(ierr);
  ierr = PetscOptionsInt("-debug","The debugging level","cimply.c",options->debug, &options->debug,NULL);CHKERRQ(ierr);
  run = options->runType;
   options->runType = (RunType) run;

   ierr = PetscOptionsBool("-interpolate","Generate intermediate mesh elements", "cimply.c",options->interpolate,&options->interpolate,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsBool("-simplex","use simplex or tensor product cells","cimply.c",options->simplex,&options->simplex,NULL);CHKERRQ(ierr);
       ierr = PetscOptionsBool("-verbose","Output additional information.","cimply.c",options->verbose,&options->verbose,NULL);CHKERRQ(ierr);
        ierr = PetscOptionsEnd(); CHKERRQ(ierr);
  return(0);
    }


 int main(int argc, char *argv[])
{
  SNES snes;			/* nonlinear solver */
  KSP ksp;                        /* the linear sovler context */
  DM dm, distributeddm;			/* problem definition */
  Vec u,r;			/* solution and residual vectors */
  Mat A,J,P;			/* Jacobian Matrix */
  AppCtx user;			/* user-defined work context */
  PetscViewer viewer;
  TS ts;                          /* in case of transient analysis */
  int dim;                        /* dimension of the anlysis */
  PetscErrorCode ierr;
  PetscMPIInt rank, numProcs;
  PetscPartitioner part;

  ierr = PetscInitialize(&argc,&argv,NULL,NULL);
  ierr = ProcessOptions(PETSC_COMM_WORLD,&user);
  ierr = DMPlexCreateFromFile(PETSC_COMM_WORLD,"nakamura.msh", PETSC_TRUE,&(dm));
  ierr = DMGetDimension(dm,&(dim));

  /* distribute the mesh */
  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
  MPI_Comm_size(PETSC_COMM_WORLD, &numProcs);
  DMPlexGetPartitioner(dm, &part);
  PetscPartitionerSetType(part, PETSCPARTITIONERPARMETIS);

  ierr = DMPlexDistribute(dm,0,NULL,&distributeddm); 
  if (distributeddm) {
    ierr=DMDestroy(&(dm));
    dm = distributeddm;
    ierr = DMDestroy(&distributeddm);
  } 

  ierr = DMView(dm,PETSC_VIEWER_STDOUT_WORLD);

  PetscPrintf(PETSC_COMM_WORLD,"starting transient analysis. Total time: %g s \n",user.time.total);
  ierr = TSCreate(PETSC_COMM_WORLD, &ts); CHKERRQ(ierr);
  ierr = TSSetType(ts, TSBEULER); CHKERRQ(ierr);
  ierr = TSSetDM(ts,dm); CHKERRQ(ierr);
  ierr = DMSetApplicationContext(dm, &user);CHKERRQ(ierr);
  ierr = SetupDiscretization(dm,&user);CHKERRQ(ierr);

  ierr = DMPlexCreateClosureIndex(dm,NULL);CHKERRQ(ierr);
  
  ierr = DMSetVecType(dm,VECMPI);CHKERRQ(ierr);    
  ierr = DMCreateGlobalVector(dm, &(u)); CHKERRQ(ierr);
  ierr = DMSetMatType(dm, MATAIJ);CHKERRQ(ierr);

  ierr = TSMonitorSet(ts,Monitor,&user,NULL);CHKERRQ(ierr);
  ierr = DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &user); CHKERRQ(ierr);
  ierr = DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &user); CHKERRQ(ierr);
  ierr = DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &user); CHKERRQ(ierr);
  

				     

  ierr = TSSetSolution(ts,u); CHKERRQ(ierr);
  ierr = TSSetDuration(ts,1000,user.time.total); CHKERRQ(ierr);
  ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP); CHKERRQ(ierr);
  ierr = TSSetEquationType(ts,TS_EQ_IMPLICIT);CHKERRQ(ierr);
  ierr = TSSetInitialTimeStep(ts,0.0, user.time.dt); CHKERRQ(ierr);
  ierr = TSSetFromOptions(ts); CHKERRQ(ierr);
  ierr = TSGetSNES(ts, &(snes)); CHKERRQ(ierr);
  ierr = SNESGetKSP(snes, &(ksp)); CHKERRQ(ierr);

  printf("solving ... .\n");
  ierr = TSSolve(ts,u); CHKERRQ(ierr);

  
  MatDestroy(&(J));
  VecDestroy(&(u));
  VecDestroy(&(r));
  TSDestroy(&(ts));
  DMDestroy(&(dm));
  PetscViewerDestroy(&(user.viewer));

  PetscFinalize();

  return 0;
}