C++ TSGetDM函数代码示例

static PetscErrorCode FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
{
  User              user = (User)ptr;
  PetscErrorCode    ierr;
  PetscScalar       **f;
  const PetscScalar **x;
  DM                dm;
  PetscInt          i,xs,xm;

  PetscFunctionBeginUser;
  if (user->reactions) {
    ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);
    ierr = DMDAVecGetArrayDOFRead(dm,X,&x);CHKERRQ(ierr);
    ierr = DMDAVecGetArrayDOF(dm,F,&f);CHKERRQ(ierr);
    ierr = DMDAGetCorners(dm,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr);

    for (i=xs; i<xs+xm; i++) {
      ierr = PetscMemcpy(user->tchemwork,x[i],(user->Nspec+1)*sizeof(x[xs][0]));CHKERRQ(ierr);
      user->tchemwork[0] *= user->Tini; /* Dimensionalize */
      ierr = TC_getSrc(user->tchemwork,user->Nspec+1,f[i]);TCCHKERRQ(ierr);
      f[i][0] /= user->Tini;           /* Non-dimensionalize */
    }

    ierr = DMDAVecRestoreArrayDOFRead(dm,X,&x);CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayDOF(dm,F,&f);CHKERRQ(ierr);
  } else {
    ierr = VecZeroEntries(F);CHKERRQ(ierr);
  }
  if (user->diffusion) {
    ierr = FormDiffusionFunction(ts,t,X,F,ptr);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

/*
   Solution - Computes the exact solution at a given time.

   Input Parameters:
   t - current time
   solution - vector in which exact solution will be computed
   appctx - user-defined application context

   Output Parameter:
   solution - vector with the newly computed exact solution
*/
PetscErrorCode Solution(TS ts,PetscReal t,Vec U,AppCtx *appctx)
{
  PetscScalar    *u,ex1,ex2,sc1,sc2,h;
  PetscErrorCode ierr;
  PetscInt       i,mstart,mend,xm,M;
  DM             da;

  ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
  ierr = DMDAGetCorners(da,&mstart,0,0,&xm,0,0);CHKERRQ(ierr);
  ierr = DMDAGetInfo(da,PETSC_IGNORE,&M,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
  h    = 1.0/M;
  mend = mstart + xm;
  /*
     Get a pointer to vector data.
  */
  ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr);

  /*
     Simply write the solution directly into the array locations.
     Alternatively, we culd use VecSetValues() or VecSetValuesLocal().
  */
  ex1 = PetscExpScalar(-36.*PETSC_PI*PETSC_PI*appctx->d*t);
  ex2 = PetscExpScalar(-4.*PETSC_PI*PETSC_PI*appctx->d*t);
  sc1 = PETSC_PI*6.*h;                 sc2 = PETSC_PI*2.*h;
  for (i=mstart; i<mend; i++) u[i] = PetscSinScalar(sc1*(PetscReal)i + appctx->a*PETSC_PI*6.*t)*ex1 + 3.*PetscSinScalar(sc2*(PetscReal)i + appctx->a*PETSC_PI*2.*t)*ex2;

  /*
     Restore vector
  */
  ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr);
  return 0;
}

PetscErrorCode FormInitialSolution(TS ts,Vec X,void *ctx)
{
  PetscScalar    **x,*xc;
  PetscErrorCode ierr;
  struct {const char *name; PetscReal massfrac;} initial[] = {
    {"CH4", 0.0948178320887},
    {"O2", 0.189635664177},
    {"N2", 0.706766236705},
    {"AR", 0.00878026702874}
  };
  PetscInt       i,j,xs,xm;
  DM             dm;

  PetscFunctionBeginUser;
  ierr = VecZeroEntries(X);CHKERRQ(ierr);
  ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);
  ierr = DMDAGetCorners(dm,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr);

  ierr = DMDAGetCoordinateArray(dm,&xc);CHKERRQ(ierr);
  ierr = DMDAVecGetArrayDOF(dm,X,&x);CHKERRQ(ierr);
  for (i=xs; i<xs+xm; i++) {
    x[i][0] = 1.0 + .05*PetscSinScalar(2.*PETSC_PI*xc[i]);  /* Non-dimensionalized by user->Tini */
    for (j=0; j<sizeof(initial)/sizeof(initial[0]); j++) {
      int ispec = TC_getSpos(initial[j].name, strlen(initial[j].name));
      if (ispec < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Could not find species %s",initial[j].name);
      ierr = PetscPrintf(PETSC_COMM_SELF,"Species %d: %s %g\n",j,initial[j].name,initial[j].massfrac);CHKERRQ(ierr);
      x[i][1+ispec] = initial[j].massfrac;
    }
  }
  ierr = DMDAVecRestoreArrayDOF(dm,X,&x);CHKERRQ(ierr);
  ierr = DMDARestoreCoordinateArray(dm,&xc);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*
   InitialConditions - Computes the solution at the initial time.

   Input Parameter:
   u - uninitialized solution vector (global)
   appctx - user-defined application context

   Output Parameter:
   u - vector with solution at initial time (global)
*/
PetscErrorCode InitialConditions(TS ts,Vec U,AppCtx *appctx)
{
  PetscScalar    *u;
  PetscErrorCode ierr;
  PetscInt       i,mstart,mend,um,M;
  DM             da;
  PetscReal      h;

  ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
  ierr = DMDAGetCorners(da,&mstart,0,0,&um,0,0);CHKERRQ(ierr);
  ierr = DMDAGetInfo(da,PETSC_IGNORE,&M,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
  h    = 1.0/M;
  mend = mstart + um;
  /*
    Get a pointer to vector data.
    - For default PETSc vectors, VecGetArray() returns a pointer to
      the data array.  Otherwise, the routine is implementation dependent.
    - You MUST call VecRestoreArray() when you no longer need access to
      the array.
    - Note that the Fortran interface to VecGetArray() differs from the
      C version.  See the users manual for details.
  */
  ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr);

  /*
     We initialize the solution array by simply writing the solution
     directly into the array locations.  Alternatively, we could use
     VecSetValues() or VecSetValuesLocal().
  */
  for (i=mstart; i<mend; i++) u[i] = PetscSinReal(PETSC_PI*i*6.*h) + 3.*PetscSinReal(PETSC_PI*i*2.*h);

  /* Restore vector */
  ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr);
  return 0;
}

/*
   Solution - Computes the exact solution at a given time

   Input Parameters:
   t - current time
   solution - vector in which exact solution will be computed
   appctx - user-defined application context

   Output Parameter:
   solution - vector with the newly computed exact solution
              u(x,t) = sin(6*PI*(x - a*t)) + 3 * sin(2*PI*(x - a*t))
*/
PetscErrorCode Solution(TS ts,PetscReal t,Vec U,AppCtx *appctx)
{
  PetscScalar    *u;
  PetscReal      a=appctx->a,h,PI6,PI2;
  PetscErrorCode ierr;
  PetscInt       i,mstart,mend,um,M;
  DM             da;

  ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
  ierr = DMDAGetCorners(da,&mstart,0,0,&um,0,0);CHKERRQ(ierr);
  ierr = DMDAGetInfo(da,PETSC_IGNORE,&M,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
  h    = 1.0/M;
  mend = mstart + um;

  /* Get a pointer to vector data. */
  ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr);

  /* u[i] = sin(6*PI*(x[i] - a*t)) + 3 * sin(2*PI*(x[i] - a*t)) */
  PI6 = PETSC_PI*6.;
  PI2 = PETSC_PI*2.;
  for (i=mstart; i<mend; i++) {
    u[i] = PetscSinReal(PI6*(i*h - a*t)) + 3.*PetscSinReal(PI2*(i*h - a*t));
  }

  /* Restore vector */
  ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr);
  return 0;
}

static PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ptr)
{
  User           user = (User)ptr;
  DM             da;
  DMDALocalInfo  info;
  PetscInt       i;
  Field          *x,*xdot,*f;
  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
  ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);

  /* Get pointers to vector data */
  ierr = DMDAVecGetArray(da,X,&x);CHKERRQ(ierr);
  ierr = DMDAVecGetArray(da,Xdot,&xdot);CHKERRQ(ierr);
  ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr);

  /* Compute function over the locally owned part of the grid */
  for (i=info.xs; i<info.xs+info.xm; i++) {
    f[i][0] = xdot[i][0] + user->k[0]*x[i][0] - user->k[1]*x[i][1] - user->s[0];
    f[i][1] = xdot[i][1] - user->k[0]*x[i][0] + user->k[1]*x[i][1] - user->s[1];
  }

  /* Restore vectors */
  ierr = DMDAVecRestoreArray(da,X,&x);CHKERRQ(ierr);
  ierr = DMDAVecRestoreArray(da,Xdot,&xdot);CHKERRQ(ierr);
  ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode FormInitialSolution(TS ts,Vec X,void *ctx)
{
  User           user = (User)ctx;
  DM             da;
  PetscInt       i;
  DMDALocalInfo  info;
  Field          *x;
  PetscReal      hx;
  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  ierr = TSGetDM(ts,&da);
  ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);
  hx   = 1.0/(PetscReal)info.mx;

  /* Get pointers to vector data */
  ierr = DMDAVecGetArray(da,X,&x);CHKERRQ(ierr);

  /* Compute function over the locally owned part of the grid */
  for (i=info.xs; i<info.xs+info.xm; i++) {
    PetscReal r = (i+1)*hx,ik = user->k[1] != 0.0 ? 1.0/user->k[1] : 1.0;
    x[i][0] = 1 + user->s[1]*r;
    x[i][1] = user->k[0]*ik*x[i][0] + user->s[1]*ik;
  }
  ierr = DMDAVecRestoreArray(da,X,&x);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

int TSFunction_Sundials(realtype t,N_Vector y,N_Vector ydot,void *ctx)
{
  TS             ts = (TS) ctx;
  DM             dm;
  DMTS           tsdm;
  TSIFunction    ifunction;
  MPI_Comm       comm;
  TS_Sundials    *cvode = (TS_Sundials*)ts->data;
  Vec            yy     = cvode->w1,yyd = cvode->w2,yydot = cvode->ydot;
  PetscScalar    *y_data,*ydot_data;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = PetscObjectGetComm((PetscObject)ts,&comm);CHKERRQ(ierr);
  /* Make the PETSc work vectors yy and yyd point to the arrays in the SUNDIALS vectors y and ydot respectively*/
  y_data    = (PetscScalar*) N_VGetArrayPointer(y);
  ydot_data = (PetscScalar*) N_VGetArrayPointer(ydot);
  ierr      = VecPlaceArray(yy,y_data);CHKERRABORT(comm,ierr);
  ierr      = VecPlaceArray(yyd,ydot_data);CHKERRABORT(comm,ierr);

  /* Now compute the right hand side function, via IFunction unless only the more efficient RHSFunction is set */
  ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);
  ierr = DMGetDMTS(dm,&tsdm);CHKERRQ(ierr);
  ierr = DMTSGetIFunction(dm,&ifunction,NULL);CHKERRQ(ierr);
  if (!ifunction) {
    ierr = TSComputeRHSFunction(ts,t,yy,yyd);CHKERRQ(ierr);
  } else {                      /* If rhsfunction is also set, this computes both parts and shifts them to the right */
    ierr = VecZeroEntries(yydot);CHKERRQ(ierr);
    ierr = TSComputeIFunction(ts,t,yy,yydot,yyd,PETSC_FALSE);CHKERRABORT(comm,ierr);
    ierr = VecScale(yyd,-1.);CHKERRQ(ierr);
  }
  ierr = VecResetArray(yy);CHKERRABORT(comm,ierr);
  ierr = VecResetArray(yyd);CHKERRABORT(comm,ierr);
  PetscFunctionReturn(0);
}

/* ------------------------------------------------------------------- */
PetscErrorCode
FormInitialSolution(TS ts, Vec U, void * /*ptr*/)
{
  DM da;
  PetscReal c = -30.0;
  PetscErrorCode ierr;
  PetscInt i, j, xs, ys, xm, ym, Mx, My;
  PetscScalar ** u;
  PetscReal hx, hy, x, y, r;

  PetscFunctionBeginUser;
  ierr = TSGetDM(ts, &da);
  CHKERRQ(ierr);
  ierr = DMDAGetInfo(da,
                     PETSC_IGNORE,
                     &Mx,
                     &My,
                     PETSC_IGNORE,
                     PETSC_IGNORE,
                     PETSC_IGNORE,
                     PETSC_IGNORE,
                     PETSC_IGNORE,
                     PETSC_IGNORE,
                     PETSC_IGNORE,
                     PETSC_IGNORE,
                     PETSC_IGNORE,
                     PETSC_IGNORE);
  CHKERRQ(ierr);

  hx = 1.0 / (PetscReal)(Mx - 1);
  hy = 1.0 / (PetscReal)(My - 1);

  /* Get pointers to vector data */
  ierr = DMDAVecGetArray(da, U, &u);
  CHKERRQ(ierr);

  /* Get local grid boundaries */
  ierr = DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
  CHKERRQ(ierr);

  /* Compute function over the locally owned part of the grid */
  for (j = ys; j < ys + ym; j++)
  {
    y = j * hy;
    for (i = xs; i < xs + xm; i++)
    {
      x = i * hx;
      r = PetscSqrtReal((x - .5) * (x - .5) + (y - .5) * (y - .5));
      if (r < .125)
        u[j][i] = PetscExpReal(c * r * r * r);
      else
        u[j][i] = 0.0;
    }
  }

  /* Restore vectors */
  ierr = DMDAVecRestoreArray(da, U, &u);
  CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

static PetscErrorCode TSComputeRHSJacobian_DMDA(TS ts,PetscReal ptime,Vec X,Mat *A,Mat *B,MatStructure *mstr,void *ctx)
{
  PetscErrorCode ierr;
  DM             dm;
  DMTS_DA        *dmdats = (DMTS_DA*)ctx;
  DMDALocalInfo  info;
  Vec            Xloc;
  void           *x;

  PetscFunctionBegin;
  if (!dmdats->rhsfunctionlocal) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"Corrupt context");
  ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);

  if (dmdats->rhsjacobianlocal) {
    ierr = DMGetLocalVector(dm,&Xloc);CHKERRQ(ierr);
    ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
    ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
    ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr);
    ierr = DMDAVecGetArray(dm,Xloc,&x);CHKERRQ(ierr);
    CHKMEMQ;
    ierr = (*dmdats->rhsjacobianlocal)(&info,ptime,x,*A,*B,mstr,dmdats->rhsjacobianlocalctx);CHKERRQ(ierr);
    CHKMEMQ;
    ierr = DMDAVecRestoreArray(dm,Xloc,&x);CHKERRQ(ierr);
    ierr = DMRestoreLocalVector(dm,&Xloc);CHKERRQ(ierr);
  } else SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_PLIB,"TSComputeRHSJacobian_DMDA() called without calling DMDATSSetRHSJacobian()");
  /* This will be redundant if the user called both, but it's too common to forget. */
  if (*A != *B) {
    ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

PetscErrorCode FormInitialSolution(TS ts,Vec X,void *ctx)
{
  User           user = (User)ctx;
  DM             da;
  PetscInt       i;
  DMDALocalInfo  info;
  Field          *x;
  PetscReal      hx;
  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  ierr = TSGetDM(ts,&da);
  ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);
  hx   = 1.0/(PetscReal)(info.mx-1);

  /* Get pointers to vector data */
  ierr = DMDAVecGetArray(da,X,&x);CHKERRQ(ierr);

  /* Compute function over the locally owned part of the grid */
  for (i=info.xs; i<info.xs+info.xm; i++) {
    PetscReal xi = i*hx;
    x[i].u = user->uleft*(1.-xi) + user->uright*xi + sin(2.*PETSC_PI*xi);
    x[i].v = user->vleft*(1.-xi) + user->vright*xi;
  }
  ierr = DMDAVecRestoreArray(da,X,&x);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

static PetscErrorCode FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
{
  User           user = (User)ptr;
  DM             da;
  DMDALocalInfo  info;
  PetscInt       i;
  PetscReal      hx;
  Field          *x,*f;
  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
  ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);
  hx   = 1.0/(PetscReal)(info.mx-1);

  /* Get pointers to vector data */
  ierr = DMDAVecGetArray(da,X,&x);CHKERRQ(ierr);
  ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr);

  /* Compute function over the locally owned part of the grid */
  for (i=info.xs; i<info.xs+info.xm; i++) {
    PetscScalar u = x[i].u,v = x[i].v;
    f[i].u = hx*(user->A + u*u*v - (user->B+1)*u);
    f[i].v = hx*(user->B*u - u*u*v);
  }

  /* Restore vectors */
  ierr = DMDAVecRestoreArray(da,X,&x);CHKERRQ(ierr);
  ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*
   Applies the second order centered difference diffusion operator on a one dimensional periodic domain
*/
static PetscErrorCode FormDiffusionFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
{
  User              user = (User)ptr;
  PetscErrorCode    ierr;
  PetscScalar       **f;
  const PetscScalar **x;
  DM                dm;
  PetscInt          i,xs,xm,j,dof;
  Vec               Xlocal;
  PetscReal         idx;

  PetscFunctionBeginUser;
  ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);
  ierr = DMDAGetInfo(dm,NULL,NULL,NULL,NULL,NULL,NULL,NULL,&dof,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
  ierr = DMGetLocalVector(dm,&Xlocal);CHKERRQ(ierr);
  ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xlocal);CHKERRQ(ierr);
  ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xlocal);CHKERRQ(ierr);
  ierr = DMDAVecGetArrayDOFRead(dm,Xlocal,&x);CHKERRQ(ierr);
  ierr = DMDAVecGetArrayDOF(dm,F,&f);CHKERRQ(ierr);
  ierr = DMDAGetCorners(dm,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr);

  idx = 1.0*user->diffus/user->dx;
  for (i=xs; i<xs+xm; i++) {
    for (j=0; j<dof; j++) {
      f[i][j] += idx*(x[i+1][j] - 2.0*x[i][j] + x[i-1][j]);
    }
  }
  ierr = DMDAVecRestoreArrayDOFRead(dm,Xlocal,&x);CHKERRQ(ierr);
  ierr = DMDAVecRestoreArrayDOF(dm,F,&f);CHKERRQ(ierr);
  ierr = DMRestoreLocalVector(dm,&Xlocal);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*
   Produces the second order centered difference diffusion operator on a one dimensional periodic domain
*/
static PetscErrorCode FormDiffusionJacobian(TS ts,PetscReal t,Vec X,Mat Amat,Mat Pmat,void *ptr)
{
  User              user = (User)ptr;
  PetscErrorCode    ierr;
  DM                dm;
  PetscInt          i,xs,xm,j,dof;
  PetscReal         idx,values[3];
  MatStencil        row,col[3];

  PetscFunctionBeginUser;
  ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);
  ierr = DMDAGetInfo(dm,NULL,NULL,NULL,NULL,NULL,NULL,NULL,&dof,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
  ierr = DMDAGetCorners(dm,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr);

  idx = 1.0*user->diffus/user->dx;
  values[0] = idx;
  values[1] = -2.0*idx;
  values[2] = idx;
  for (i=xs; i<xs+xm; i++) {
    for (j=0; j<dof; j++) {
      row.i = i;      row.c = j;
      col[0].i = i-1; col[0].c = j;
      col[1].i = i;   col[1].c = j;
      col[2].i = i+1; col[2].c = j;
      ierr = MatSetValuesStencil(Pmat,1,&row,3,col,values,ADD_VALUES);CHKERRQ(ierr);
    }
  }
  ierr = MatAssemblyBegin(Pmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(Pmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*
 * This is a modified version of PETSc/src/ts/examples/tutorials/ex15.c
 * to demonstrate how MOOSE interact with an external solver package
 */
PetscErrorCode
externalPETScDiffusionFDMSolve(TS ts, Vec u, PetscReal dt, PetscReal time)
{
  PetscErrorCode ierr;
#if !PETSC_VERSION_LESS_THAN(3, 8, 0)
  PetscInt current_step;
#endif
  DM da;

  PetscFunctionBeginUser;

  ierr = TSGetDM(ts, &da);
  CHKERRQ(ierr);

#if !PETSC_VERSION_LESS_THAN(3, 7, 0)
  PetscOptionsSetValue(NULL, "-ts_monitor", NULL);
  PetscOptionsSetValue(NULL, "-snes_monitor", NULL);
  PetscOptionsSetValue(NULL, "-ksp_monitor", NULL);
#else
  PetscOptionsSetValue("-ts_monitor", NULL);
  PetscOptionsSetValue("-snes_monitor", NULL);
  PetscOptionsSetValue("-ksp_monitor", NULL);
#endif

  /*ierr = TSSetMaxTime(ts,1.0);CHKERRQ(ierr);*/
  ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
  CHKERRQ(ierr);

  ierr = TSSetSolution(ts, u);
  CHKERRQ(ierr);
  ierr = TSSetTimeStep(ts, dt);
  CHKERRQ(ierr);
  ierr = TSSetTime(ts, time - dt);
  CHKERRQ(ierr);
#if !PETSC_VERSION_LESS_THAN(3, 8, 0)
  ierr = TSGetStepNumber(ts, &current_step);
  CHKERRQ(ierr);
  ierr = TSSetMaxSteps(ts, current_step + 1);
  CHKERRQ(ierr);
#else
  SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Require PETSc-3.8.x or higher ");
#endif
  /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
   Sets various TS parameters from user options
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  ierr = TSSetFromOptions(ts);
  CHKERRQ(ierr);
  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Solve nonlinear system
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  ierr = TSSolve(ts, u);
  CHKERRQ(ierr);
  PetscFunctionReturn(0);
}