#ifdef PETSC_RCS_HEADER
static char vcid[] = "$Id: gts.c,v 1.9 2000/01/10 03:54:19 knepley Exp $";
#endif

#include "src/ts/tsimpl.h"        /*I "petscts.h" I*/
#include "src/gvec/gvecimpl.h"    /*I "gsolver.h" I*/
#include "gsolver.h"

/* Logging support */
int GTS_Reform, GTS_Reallocate;

/*
  This file provides routines for grid TS objects. They support solution of
  systems of time-dependent nonlinear equations generated from a mesh and
  discretization specified by a Grid object.
*/

#undef  __FUNCT__
#define __FUNCT__ "GTSDestroy"
/*@C
  GTSDestroy - Destroys a grid TS.

  Collective on GTS

  Input Parameters:
. ts - The nonlinear solver context

  Level: beginner

.keywords: grid TS, destroy
.seealso: TSDestroy(), GTSDuplicate()
@*/
int GTSDestroy(GTS ts)
{
  TSProblemType type;
  int           ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ts, TS_COOKIE);
  if (--ts->refct > 0)
    PetscFunctionReturn(0);
  ierr = TSGetProblemType(ts, &type);                                                                     CHKERRQ(ierr);
  if (type == TS_LINEAR) {
    ierr = GMatDestroy(ts->A);                                                                            CHKERRQ(ierr);
  }
  ierr = PetscFree(ts->isExplicit);                                                                       CHKERRQ(ierr);
  ierr = PetscFree(ts->Iindex);                                                                           CHKERRQ(ierr);
  if (ts->snes) {
    ierr = GSNESDestroy(ts->snes);                                                                        CHKERRQ(ierr);
    ts->snes = PETSC_NULL;
  }
  ierr = TSDestroy(ts);                                                                                   CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GTSView"
/*@ 
  GTSView - Views a grid TS.

  Collective on GTS

  Input Parameters:
. ts     - The grid TS
. viewer - An optional visualization context

  Options Database Key:
$ -ts_view : calls TSView() at end of TSSolve()

  Notes:
  The available visualization contexts include
$   VIEWER_STDOUT_SELF - standard output (default)
$   VIEWER_STDOUT_WORLD - synchronized standard
$     output where only the first processor opens
$     the file.  All other processors send their 
$     data to the first processor to print. 

  The user can open alternative vistualization contexts with
$   PetscViewerFileOpenASCII() - output vector to a specified file

  Level: beginner

.keywords: view, visualize, output, print, write, draw
.seealso: TSView()
@*/
int GTSView(GTS ts, PetscViewer viewer)
{  
  Grid       grid;
  FILE      *fd;
  int        numActiveFields;
  int        field, func, op;
  PetscTruth isascii;
  int        ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ts, TS_COOKIE);
  if (!viewer) {
    viewer = PETSC_VIEWER_STDOUT_SELF;
  } else {
    PetscValidHeaderSpecific(viewer, PETSC_VIEWER_COOKIE);
  }

  ierr = GTSGetGrid(ts, &grid);                                                                           CHKERRQ(ierr);
  ierr = PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_ASCII, &isascii);                            CHKERRQ(ierr);
  if (isascii == PETSC_TRUE) {
    ierr = GridView(grid, viewer);                                                                        CHKERRQ(ierr);
    ierr = PetscViewerFlush(viewer);                                                                      CHKERRQ(ierr);
    ierr = PetscViewerASCIIGetPointer(viewer, &fd);                                                       CHKERRQ(ierr);
    PetscFPrintf(ts->comm, fd, "GTS Object:\n");
    ierr = GridGetNumActiveFields(grid, &numActiveFields);                                                CHKERRQ(ierr);
    if (numActiveFields == 1) {
      PetscFPrintf(ts->comm, fd, "  %d active field:\n", numActiveFields);
    } else {
      PetscFPrintf(ts->comm, fd, "  %d active fields:\n", numActiveFields);
    }
    for(field = 0; field < grid->numFields; field++) {
      if (grid->fields[field].isActive == PETSC_FALSE) continue;
      if (grid->fields[field].name) {
        PetscFPrintf(grid->comm, fd, "  %s field with %d components:\n", grid->fields[field].name, grid->fields[field].numComp);
      } else {
        PetscFPrintf(grid->comm, fd, "  field %d with %d components:\n", field, grid->fields[field].numComp);
      }
      if (ts->isExplicit[field]) PetscFPrintf(ts->comm, fd, "    Explicitly time dependent\n");
      for(func = 0; func < grid->numRhsFuncs; func++) {
        if (grid->rhsFuncs[func].field == field) {
          PetscFPrintf(grid->comm, fd, "    Rhs function with scalar multiplier %g\n", PetscRealPart(grid->rhsFuncs[func].alpha));
        }
      }
      for(op = 0; op < grid->numRhsOps; op++) {
        if (grid->rhsOps[op].field == field) {
          if (grid->rhsOps[op].nonlinearOp != PETSC_NULL) {
            PetscFPrintf(grid->comm, fd, "    Rhs nonlinear operator with scalar multiplier %g\n",
                         PetscRealPart(grid->rhsOps[field].alpha));
          } else {
            PetscFPrintf(grid->comm, fd, "    Rhs operator %d with scalar multiplier %g\n",
                         grid->rhsOps[op].op, PetscRealPart(grid->rhsOps[op].alpha));
          }
        }
      }
      for(op = 0; op < grid->numMatOps; op++) {
        if (grid->matOps[op].field == field) {
          PetscFPrintf(grid->comm, fd, "    Jacobian operator %d with scalar multiplier %g\n",
                       grid->matOps[op].op, PetscRealPart(grid->matOps[op].alpha));
        }
      }
    }
    ierr = TSView(ts, viewer);                                                                            CHKERRQ(ierr);
  }

  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GTSSerialize"
/*@ 
  GTSSerialize - This function stores or recreates a grid TS using a viewer for a binary file.

  Collective on grid

  Input Parameters:
. comm   - The communicator for the grid TS object
. viewer - The viewer context
. store  - This flag is PETSC_TRUE is data is being written, otherwise it will be read

  Output Parameter:
. ts     - The grid TS

  Level: beginner

.keywords: grid TS, serialize
.seealso: GridSerialize()
@*/
int GTSSerialize(Grid grid, GTS *ts, PetscViewer viewer, PetscTruth store)
{
  MPI_Comm comm;
  int      ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid,   GRID_COOKIE);
  PetscValidHeaderSpecific(viewer, PETSC_VIEWER_COOKIE);
  PetscValidPointer(ts);
#ifndef PETSC_USE_DYNAMIC_LIBRARIES
  ierr = GSolverInitializePackage(PETSC_NULL);                                                            CHKERRQ(ierr);
#endif


  ierr = PetscObjectGetComm((PetscObject) grid, &comm);                                                   CHKERRQ(ierr);
  /* We need a way to communicate the grid to the TS */
  if (store == PETSC_FALSE) {
    *ts = (GTS) grid;
  }
  ierr = TSSerialize(comm, ts, viewer, store);                                                            CHKERRQ(ierr);
  /* Newly created object should have grid as a parent */
  if (store == PETSC_FALSE) {
    ierr = PetscObjectCompose((PetscObject) *ts, "Grid", (PetscObject) grid);                             CHKERRQ(ierr);

    /* Set boundary conditions */
    ierr = TSSetRhsBC(*ts, GTSRhsBC);                                                                     CHKERRQ(ierr);
    ierr = TSSetSolutionBC(*ts, GTSSolutionBC);                                                           CHKERRQ(ierr);
    /* Save space/time by reducing the system at the element level */
    if ((*ts)->problem_type == TS_NONLINEAR) {
      ierr = GridSetReduceElement(grid, PETSC_TRUE);                                                      CHKERRQ(ierr);
    }

    /* Set update functions */
    ierr = TSSetPreStep(*ts, GTSPreStep);                                                                 CHKERRQ(ierr);
    ierr = TSSetUpdate(*ts, GTSUpdate);                                                                   CHKERRQ(ierr);
    ierr = TSSetPostStep(*ts, GTSPostStep);                                                               CHKERRQ(ierr);
  }

  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GTSDuplicate"
/*@C
  GTSDuplicate - Duplicates a grid TS.

  Collective on GTS

  Input Parameter:
. ts    - The grid TS

  Output Parameter:
. newts - The new grid TS

  Level: beginner

.keywords: grid ts, duplicate
.seealso: TSDuplicate(), GTSDestroy()
@*/
int GTSDuplicate(GTS ts, GTS *newts)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ts, TS_COOKIE);
  SETERRQ(PETSC_ERR_SUP, " ");
#if 0
  PetscFunctionReturn(0);
#endif
}

#undef  __FUNCT__
#define __FUNCT__ "GTSReform"
/*@C
  GTSReform - Cleans up any auxilliary structures in the GTS. Should
  be called after GridReform(), but before any user specific changes
  to the grid. Usually followed by GTSReallocate().

  Collective on GTS

  Input Parameter:
. ts - The grid TS

  Level: advanced

.keywords: grid ts, reform
.seealso: GTSReallocate(), GridReform()
@*/
int GTSReform(GTS ts)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ts, TS_COOKIE);
  if (ts->ops->reform) {
    ierr = PetscLogEventBegin(GTS_Reform, ts, 0, 0, 0);                                                     CHKERRQ(ierr);
    ierr = (*ts->ops->reform)(ts);                                                                        CHKERRQ(ierr);
    ierr = PetscLogEventEnd(GTS_Reform, ts, 0, 0, 0);                                                       CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GTSReallocate"
/*@C
  GTSReallocate - Reallocates the storage for a grid TS. Should be
  called after GTSReform() and any user updating of the grid.

  Collective on GTS

  Input Parameters:
. ts - The grid TS
. x  - The new solution vector

  Level: advanced

.keywords: grid ts, reallocate
.seealso: GTSReform(), GridReform()
@*/
int GTSReallocate(GTS ts, GVec x)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ts, TS_COOKIE);
  PetscValidHeaderSpecific(x,  VEC_COOKIE);
  ierr = TSSetSolution(ts, x);                                                                            CHKERRQ(ierr);
  if (ts->ops->reallocate) {
    ierr = PetscLogEventBegin(GTS_Reallocate, ts, 0, 0, 0);                                               CHKERRQ(ierr);
    ierr = (*ts->ops->reallocate)(ts);                                                                    CHKERRQ(ierr);
    ierr = PetscLogEventEnd(GTS_Reallocate, ts, 0, 0, 0);                                                 CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GTSGetGrid"
/*@
  GTSGetGrid - Gets the grid from a grid TS.

  Not collective

  Input Parameter:
. ts   - The grid TS

  Output Parameter:
. grid - The grid context

  Level: intermediate

.keywords: grid ts, get
.seealso: GVecGetGrid(), GMatGetGrid(), GSNESGetGrid()
@*/
int GTSGetGrid(GTS ts, Grid *grid)
{
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ts, TS_COOKIE);
  PetscValidPointer(grid);

  ierr = PetscObjectQuery((PetscObject) ts, "Grid", (PetscObject *) grid);                                CHKERRQ(ierr);
  PetscValidHeaderSpecific(*grid, GRID_COOKIE);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GTSGetInitialTimeStep"
/*@
  GTSGetInitialTimeStep - Returns the initial time step

  Not collective

  Input Parameter:
. ts - The grid TS

  Output Parameter:
. dt - The initial time step

  Level: intermediate

.keywords: grid ts, get, time step
.seealso: TSSetInitialTimeStep()
@*/
int GTSGetInitialTimeStep(GTS ts, double *dt)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ts, TS_COOKIE);
  PetscValidScalarPointer(dt);
  *dt = ts->initial_time_step;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GTSGetTimeDependence"
/*@
  GTSGetTimeDependence - Determines whether or not a field has
  explicit time dependence.

  Not collective

  Input Parameters:
. ts    - The grid TS
. field - The field

  Output Parameter:
. flag  - Indicates whether or not the field has explicit time dependence

  Level: intermediate

.keywords: grid ts, get
.seealso: GTSSetTimeDependence()
@*/
int GTSGetTimeDependence(GTS ts, int field, PetscTruth *flag)
{
  Grid grid;
  int  numFields;
  int  ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ts, TS_COOKIE);
  PetscValidPointer(flag);
  ierr = GTSGetGrid(ts, &grid);                                                                           CHKERRQ(ierr);
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  ierr = GridGetNumFields(grid, &numFields);                                                              CHKERRQ(ierr);
  if ((field < 0) || (field > numFields)) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Invalid field number");
  *flag = ts->isExplicit[field];
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GTSSetContext"
/*@
  GTSSetContext - Sets a user-defined context for use in calculation.

  Collective on GTS

  Input Parameters:
. ts  - The grid TS
. ctx - The context

  Level: intermediate

.keywords: grid ts, context
.seealso: GTSGetContext()
@*/
int GTSSetContext(GTS ts, void *ctx)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ts, TS_COOKIE);
  ts->funP = ctx;
  ts->jacP = ctx;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GTSSetTimeDependence"
/*@
  GTSSetTimeDependence - Sets the explicit time dependence of a field.

  Collective on GTS

  Input Parameters:
. ts    - The grid TS
. field - The field
. flag  - Indicates whether or not the field has explicit time dependence

  Level: intermediate

.keywords: grid ts, get
.seealso: GTSGetTimeDependence()
@*/
int GTSSetTimeDependence(GTS ts, int field, PetscTruth flag)
{
  Grid grid;
  int  numFields;
  int  ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ts, TS_COOKIE);
  ierr = GTSGetGrid(ts, &grid);                                                                           CHKERRQ(ierr);
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  ierr = GridGetNumFields(grid, &numFields);                                                              CHKERRQ(ierr);
  if ((field < 0) || (field > numFields)) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Invalid field number");
  ts->isExplicit[field] = flag;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GTSEvaluateRhs"
/*@
  GTSEvaluateRhs - Constructs the Rhs vector for U_t = F(U, t).

  Collective on GTS

  Input Parameters:
. ts   - The grid TS
. x    - The current iterate
. ctx  - The optional user context

  Output Parameter:
. f    - The function value

  Level: advanced

.keywords: grid ts, rhs
.seealso: GSNESEvaluateRhs()
@*/
int GTSEvaluateRhs(GTS ts, double t, GVec x, GVec f, PetscObject ctx)
{
  Grid        grid;
  PetscObject oldCtx;
  PetscScalar zero  = 0.0;
  int         ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ts, TS_COOKIE);
  ierr = GTSGetGrid(ts, &grid);                                                                           CHKERRQ(ierr);
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  /* Initialize vector */
  ierr = VecSet(&zero, f);                                                                                CHKERRQ(ierr);
  /* Setup new context */
  ierr = GTSCreateContext(ts, t, ctx, &oldCtx);                                                           CHKERRQ(ierr);
  /* Form function */
  ierr = GridEvaluateRhs(grid, x, f, ctx);                                                                CHKERRQ(ierr);
  /* Cleanup */
  ierr = GTSDestroyContext(ts, ctx, oldCtx);                                                              CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GTSEvaluateJacobian"
/*@
  GTSEvaluateJacobian - Constructs the Jacobian matrix for U_t = F(U,t)

  Collective on GTS

  Input Parameters:
. ts   - The grid TS
. x    - The current iterate
. flag - The matrix structure flag
. ctx  - The optional user context

  Output Parameters:
. J    - The Jacobian matrix
. M    - The preconditioner for the Jacobian matrix, usually the same as J

  Level: advanced

.keywords: grid ts, evaluate, jacobian
.seealso: SNESComputeJacobian(), GridEvaluateSystemMatrix()
@*/
int GTSEvaluateJacobian(GTS ts, double t, GVec x, GMat *J, GMat *M, MatStructure *flag, PetscObject ctx)
{
  Grid        grid;
  PetscObject oldCtx;
  int         ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ts, TS_COOKIE);
  ierr = GTSGetGrid(ts, &grid);                                                                           CHKERRQ(ierr);
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  /* Initialize matrix */
  ierr = MatZeroEntries(*J);                                                                              CHKERRQ(ierr);
  /* Setup new context */
  ierr = GTSCreateContext(ts, t, ctx, &oldCtx);                                                           CHKERRQ(ierr);
  /* Form function */
  ierr = GridEvaluateSystemMatrix(grid, x, J, M, flag, ctx);                                              CHKERRQ(ierr);
  /* Apply boundary conditions */
  ierr = GMatSetBoundary(*J, 1.0, ctx);                                                                   CHKERRQ(ierr);
  /* Cleanup */
  ierr = GTSDestroyContext(ts, ctx, oldCtx);                                                              CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GTSEvaluateSystemMatrix"
/*@
  GTSEvaluateSystemMatrix - Constructs the system matrix for U_t = A(U,t) U

  Collective on GTS

  Input Parameters:
. ts   - The grid TS
. flag - The matrix structure flag
. ctx  - The optional user context

  Output Parameters:
. A    - The system matrix
. M    - The preconditioner for the system matrix, usually the same as A

  Level: advanced

.keywords: grid ts, evaluate, jacobian
.seealso: SNESComputeJacobian(), GridEvaluateSystemMatrix()
@*/
int GTSEvaluateSystemMatrix(GTS ts, double t, GMat *A, GMat *M, MatStructure *flag, PetscObject ctx)
{
  Grid        grid;
  PetscObject oldCtx;
  int         ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ts, TS_COOKIE);
  ierr = GTSGetGrid(ts, &grid);                                                                           CHKERRQ(ierr);
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  /* Initialize matrix */
  ierr = MatZeroEntries(*A);                                                                              CHKERRQ(ierr);
  /* Setup new context */
  ierr = GTSCreateContext(ts, t, ctx, &oldCtx);                                                           CHKERRQ(ierr);
  /* Form function */
  ierr = GridEvaluateSystemMatrix(grid, PETSC_NULL, A, M, flag, ctx);                                     CHKERRQ(ierr);
  /* Apply boundary conditions */
  ierr = GMatSetBoundary(*A, 1.0, ctx);                                                                   CHKERRQ(ierr);
  /* Cleanup */
  ierr = GTSDestroyContext(ts, ctx, oldCtx);                                                              CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GridCalcBCValues"
/*@
  GTSCalcBCValues - This function calculates the boundary values. It
  is normally called once a timestep when using time dependent boundary
  conditions.

  Collective on GTS

  Input Parameter:
. ts - The GTS

  Level: advanced

.keywords: grid ts, reduction, boundary conditions
.seealso: GridSetBCContext(), GridSetBC(), GridAddBC()
@*/
int GTSCalcBCValues(GTS ts)
{
  Grid        grid;
  PetscObject oldCtx;
  int         ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ts, TS_COOKIE);
  ierr = GTSGetGrid(ts, &grid);                                                                           CHKERRQ(ierr);
  /* Setup new context */
  ierr = GTSCreateContext(ts, ts->ptime, (PetscObject) ts->funP, &oldCtx);                                CHKERRQ(ierr);
  /* Form function */
  ierr = GridCalcBCValues(grid, PETSC_TRUE, ts->funP);                                                    CHKERRQ(ierr);
  /* Cleanup */
  ierr = GTSDestroyContext(ts, (PetscObject) ts->funP, oldCtx);                                           CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "GTSRhsBC"
/*@
  GTSRhsBC - This sets the boundary conditions registered with the grid
  on the Rhs vector, or zero boundary conditions for nonlinear systems.

  Collective on GTS

  Input Parameters:
. ts  - The TS context obtained from TSCreate()
. ctx - The user-supplied context

  Output Parameter:
. rhs - The Rhs

  Level: advanced

.keywords: grid ts, timestep, set, boundary conditions
.seealso: GridSetBCContext(), GridSetBC(), GridAddBC()
@*/
int GTSRhsBC(GTS ts, GVec rhs, void *ctx)
{
  TSProblemType type;
  PetscObject   oldCtx;
  int           ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ts,  TS_COOKIE);
  PetscValidHeaderSpecific(rhs, VEC_COOKIE);
  ierr = TSGetProblemType(ts, &type);                                                                     CHKERRQ(ierr);
  /* Setup new context */
  ierr = GTSCreateContext(ts, ts->ptime, (PetscObject) ctx, &oldCtx);                                     CHKERRQ(ierr);
  switch(type)
  {
  case TS_LINEAR:
    ierr = GVecSetBoundary(rhs, ctx);                                                                     CHKERRQ(ierr);
    break;
  case TS_NONLINEAR:
    ierr = GVecSetBoundaryZero(rhs, ctx);                                                                 CHKERRQ(ierr);
    /* ierr = TSGetSolution(ts, &u);                                                                      CHKERRQ(ierr);        */
    /* ierr = GVecSetBoundaryDifference(rhs, u, &GTSctx);                                                 CHKERRQ(ierr); */
    break;
  default:
    SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid TS problem type %d", type);
  }
  /* Cleanup */
  ierr = GTSDestroyContext(ts, (PetscObject) ctx, oldCtx);                                                CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "GTSSolutionBC"
/*@
  GTSSolutionBC - This sets the boundary conditions registered with the grid
  on the solution vector.

  Collective on GTS

  Input Parameters:
. ts  - The TS context obtained from TSCreate()
. ctx - The user-supplied context

  Output Parameter:
. sol - The solution

  Level: advanced

.keywords: grid ts, timestep, set, boundary conditions
.seealso: GridSetBCContext(), GridSetBC(), GridAddBC()
@*/
int GTSSolutionBC(GTS ts, GVec sol, void *ctx)
{
  PetscObject oldCtx;
  int         ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ts,  TS_COOKIE);
  PetscValidHeaderSpecific(sol, VEC_COOKIE);
  /* Setup new context */
  ierr = GTSCreateContext(ts, ts->ptime, (PetscObject) ctx, &oldCtx);                                     CHKERRQ(ierr);
  /* Apply solution BC */
  ierr = GVecSetBoundary(sol, ctx);                                                                       CHKERRQ(ierr);
  /* Cleanup */
  ierr = GTSDestroyContext(ts, (PetscObject) ctx, oldCtx);                                                CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "GTSSolutionBCforGSNES"
/*@
  GTSSolutionBCforGSNES - This functions allows a GSNES
  imbedded in a GTS to correctly set time dependent boundary
  conditions on the initial solution at each time step.

  Collective on GTS

  Input Parameters:
. ts  - The TS context obtained from TSCreate()
. ctx - The user-supplied context

  Output Parameter:
. sol - The solution

  Level: advanced

.keywords: grid ts, timestep, set, boundary conditions
.seealso: GridSetBCContext(), GridSetBC(), GridAddBC()
@*/
int GTSSolutionBCforGSNES(GSNES snes, GVec sol, void *ctx)
{
  GTS ts = (GTS) ctx;
  int ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ts,  TS_COOKIE);
  PetscValidHeaderSpecific(sol, VEC_COOKIE);
  ierr = GTSSolutionBC(ts, sol, ts->funP);                                                                CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "GTSPreStep"
/*@
  GTSPreStep - This is a general purpose function which
  is called once at the beginning of time stepping. It
  currently performs context wrapping and grid manipulation.

  Collective on GTS

  Input Parameters:
. ts  - The GTS context

  Level: advanced

.keywords: grid ts, timestep
.seealso: GTSUpdate(), GTSPostStep()
@*/
int GTSPreStep(GTS ts)
{
  Grid          grid;
  PetscScalar   mdt = 1.0/ts->time_step;
  PetscTruth    isTimeDependent;
  TSProblemType type;
  int           numFields;
  int           field, f;
  int           ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ts, TS_COOKIE);

  ierr = GTSGetGrid(ts, &grid);                                                                           CHKERRQ(ierr);
  ierr = TSGetProblemType(ts, &type); CHKERRQ(ierr);
  if (type == TS_NONLINEAR) {
    /* Scale Rhs so we get u_t - F = 0 */
    ierr = GridScaleRhs(grid, -1.0);                                                                      CHKERRQ(ierr);

    /* Add the identity operator and scale to get J_{sys} = - J_{F} where J_{F} is the given Jacobian of F */
    ierr = GridScaleSystemMatrix(grid, -1.0);                                                             CHKERRQ(ierr);
  }

  /* Add the identity operator and scale to get J = \Delta t^{-1} I + J_{sys} where J_{sys} is the system matrix */
  ierr = GridGetNumActiveFields(grid, &numFields);                                                        CHKERRQ(ierr);
  for(f = 0; f < numFields; f++) {
    ierr = GridGetActiveField(grid, f, &field);                                                           CHKERRQ(ierr);
    ierr = GTSGetTimeDependence(ts, field, &isTimeDependent);                                             CHKERRQ(ierr);
    if (isTimeDependent) {
      ierr = GridAddMatOperator(grid, field, field, IDENTITY, mdt, PETSC_FALSE, &ts->Iindex[field]);      CHKERRQ(ierr);
    }
  }

  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "GTSUpdate"
/*@
  GTSUpdate - This is a general purpose function which
  is called at the beginning of every time step. It
  currently calculates the mesh acceleration and moves
  the mesh.

  Collective on GTS

  Input Parameters:
. ts  - The GTS context
.  t   - The current time

  Output Parameter:
. dt  - The current time step

  Level: advanced

.keywords: grid ts, timestep
.seealso: GTSPreStep(), GTSPostStep()
@*/
int GTSUpdate(GTS ts, PetscReal t, PetscReal *dt)
{
  Grid       grid;
  Mesh       mesh;
  MeshMover  mover;
  PetscTruth isMoving;
  int        ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ts, TS_COOKIE);

  ierr = GTSGetGrid(ts, &grid);                                                                           CHKERRQ(ierr);
  ierr = GridGetMesh(grid, &mesh);                                                                        CHKERRQ(ierr);
  ierr = MeshGetMovement(mesh, &isMoving);                                                                CHKERRQ(ierr);
  if (isMoving == PETSC_TRUE) {
    ierr = MeshGetMover(mesh, &mover);                                                                    CHKERRQ(ierr);
    /* Calculate mesh acceleration */
    ierr = MeshMoverCalcNodeAccelerations(mover, PETSC_TRUE);                                             CHKERRQ(ierr);
    /* Move the mesh */
    ierr = MeshMoveMesh(mesh, ts->time_step);                                                             CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "GTSPostStep"
/*@
  GTSPostStep - This is a general purpose function which
  is called once at the end of time stepping. It currently
  performs context unwrapping and grid manipulation.

  Collective on GTS

  Input Parameters:
. ts  - The GTS context

  Level: advanced

.keywords: grid ts, timestep
.seealso: GTSPreStep(), GTSUpdate()
@*/
int GTSPostStep(GTS ts)
{
  Grid          grid;
  TSProblemType type;
  PetscTruth    isTimeDependent;
  int           numFields;
  int           field, f;
  int           ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ts, TS_COOKIE);

  ierr = GTSGetGrid(ts, &grid);                                                                           CHKERRQ(ierr);
  ierr = TSGetProblemType(ts, &type); CHKERRQ(ierr);
  if (type == TS_NONLINEAR) {
    /* Scale Rhs back to the original state */
    ierr = GridScaleRhs(grid, -1.0);                                                                      CHKERRQ(ierr);

    /* Remove the identity and scale back */
    ierr = GridScaleSystemMatrix(grid, -1.0);                                                             CHKERRQ(ierr);
  }

  ierr = GridGetNumActiveFields(grid, &numFields);                                                        CHKERRQ(ierr);
  for(f = 0; f < numFields; f++) {
    ierr = GridGetActiveField(grid, f, &field);                                                           CHKERRQ(ierr);
    ierr = GTSGetTimeDependence(ts, field, &isTimeDependent);                                             CHKERRQ(ierr);
    if (isTimeDependent) {
      ierr = GridRemoveMatOperator(grid, ts->Iindex[field]);                                              CHKERRQ(ierr);
    }
  }

  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GTSCreate"
/*@
  GTSCreate - Creates a grid TS associated with a
  particular discretization context.

  Collective on Grid

  Input Parameter:
. grid - The discretization context
. ctx  - The optional user context

  Output Parameter:
. ts   - The integrator context

  Level: beginner

  Note:
  The TS type is TS_NONLINEAR if any nonlinear operators have
  been activated on the grid.

.keywords: grid ts, create
.seealso: TSCreate(), TSSetType()
@*/
int GTSCreate(Grid grid, void *ctx, GTS *ts)
{
  MPI_Comm      comm;
  TSProblemType type;
  GTS           gts;
  GMat          A;
  int           numFields, numRhsOps;
  int           field, op;
  int           ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(grid, GRID_COOKIE);
  PetscValidPointer(ts);
#ifndef PETSC_USE_DYNAMIC_LIBRARIES
  ierr = GSolverInitializePackage(PETSC_NULL);                                                            CHKERRQ(ierr);
#endif

  /* Get TS type */
  type = TS_LINEAR;
  ierr = GridGetNumRhsOperators(grid, &numRhsOps);                                                        CHKERRQ(ierr);
  for(op = 0; op < grid->numRhsOps; op++) {
    if (grid->rhsOps[op].nonlinearOp != PETSC_NULL) {
      type = TS_NONLINEAR;
      break;
    }
  }

  /* Create the TS */
  ierr = PetscObjectGetComm((PetscObject) grid, &comm);                                                   CHKERRQ(ierr);
  ierr = TSCreate(comm, &gts);                                                                            CHKERRQ(ierr);
  ierr = TSSetProblemType(gts, type);                                                                     CHKERRQ(ierr);
  gts->isGTS = PETSC_TRUE;
  gts->bops->destroy = (int (*)(PetscObject)) GTSDestroy;
  gts->bops->view    = (int (*)(PetscObject, PetscViewer)) GTSView;

  /* Set boundary conditions */
  ierr = TSSetRhsBC(gts, GTSRhsBC);                                                                       CHKERRQ(ierr);
  ierr = TSSetSolutionBC(gts, GTSSolutionBC);                                                             CHKERRQ(ierr);
  /* Save space/time by reducing the system at the element level */
  if (type == TS_NONLINEAR) {
    ierr = GridSetReduceElement(grid, PETSC_TRUE);                                                        CHKERRQ(ierr);
  }

  /* Set update functions */
  ierr = TSSetPreStep(gts, GTSPreStep);                                                                   CHKERRQ(ierr);
  ierr = TSSetUpdate(gts, GTSUpdate);                                                                     CHKERRQ(ierr);
  ierr = TSSetPostStep(gts, GTSPostStep);                                                                 CHKERRQ(ierr);

  /* Create the array which indicates time dependent fields */
  ierr = GridGetNumFields(grid, &numFields);                                                              CHKERRQ(ierr);
  ierr = PetscMalloc(numFields * sizeof(PetscTruth), &gts->isExplicit);                                   CHKERRQ(ierr);
  PetscLogObjectMemory(gts, numFields * sizeof(PetscTruth));
  for(field = 0; field < numFields; field++) {
    gts->isExplicit[field] = PETSC_TRUE;
  }

  /* Create the array which holds the index of the identity operator for each time dependent field */
  ierr = PetscMalloc(numFields * sizeof(int), &gts->Iindex);                                              CHKERRQ(ierr);
  PetscLogObjectMemory(gts, numFields * sizeof(int));
  for(field = 0; field < numFields; field++)
    gts->Iindex[field]   = -1;

  /* Setup TS data structures */
  gts->funP = ctx;
  gts->jacP = ctx;
  if (type == TS_LINEAR) {
    PetscTruth reduceSystem;

    ierr = GMatCreate(grid, &A);                                                                          CHKERRQ(ierr);
    ierr = GridGetReduceSystem(grid, &reduceSystem);                                                      CHKERRQ(ierr);
    if (reduceSystem == PETSC_FALSE) {
      ierr = MatSetOption(A, MAT_YES_NEW_NONZERO_LOCATIONS);                                              CHKERRQ(ierr);
    }
    gts->A = A;
    gts->B = A;
    PetscLogObjectParent(grid, A);
  }

  ierr = PetscObjectCompose((PetscObject) gts, "Grid", (PetscObject) grid);                               CHKERRQ(ierr);
  *ts = gts;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GTSSolutionMonitor"
/*@
  GTSSolutionMonitor - Monitors solution at each GTS iteration.

  Collective on GTS

  Input Parameters:
. ts    - The integrator context
. step  - The current time step
. ltime - The current time
. sol   - The current solution vector
. ctx   - The viewer

  Level: intermediate

.keywords: grid ts, monitor, solution
.seealso: TSDefaultMonitor(),TSSetMonitor(),GTSResidualMonitor(),
          GTSErrorMonitor()
@*/
int GTSSolutionMonitor(GTS ts, int step, PetscReal ltime, GVec sol, void *ctx)
{
  PetscViewer viewer = (PetscViewer) ctx;
  int         ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ts,  TS_COOKIE);
  PetscValidHeaderSpecific(sol, VEC_COOKIE);
  ierr = GVecView(sol, viewer);                                                                           CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GTSErrorMonitor"
/*@
  GTSErrorMonitor - Displays the error at each iteration.

  Collective on GTS

  Input Parameters:
+ ts      - The integrator context
. step    - The current time step
. ltime   - The current time
. sol     - The current solution vector
- monCtx  - The GTSErrorMonitorCxt

  Notes: 
  The final argument to TSSetMonitor() with this routine must be
  a pointer to a GTSErrorMonitorCtx.

  Level: intermediate

.keywords: grid ts, monitor, error
.seealso: TSDefaultMonitor(),TSSetMonitor(),GTSSolutionMonitor(), GTSResidualMonitor()
@*/
int GTSErrorMonitor(GTS ts, int step, PetscReal ltime, GVec sol, void *monCtx)
{
  GTSErrorMonitorCtx *errorCtx = (GTSErrorMonitorCtx *) monCtx;
  PetscObject         oldCtx;
  PetscScalar         minusOne = -1.0;
  GVec                e;
  MPI_Comm            comm;
  FILE               *file;
  PetscReal           norm_2, norm_max, sol_norm_2;
  int                 ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ts,  TS_COOKIE);
  PetscValidHeaderSpecific(sol, VEC_COOKIE);

  ierr = GTSCreateContext(ts, ltime, (PetscObject) errorCtx->ctx, &oldCtx);                                CHKERRQ(ierr);
  ierr = (*errorCtx->solutionFunc)(errorCtx->solution, errorCtx->ctx);                                    CHKERRQ(ierr);
  ierr = GVecGetWorkGVec(errorCtx->solution, &e);                                                         CHKERRQ(ierr);
  ierr = VecWAXPY(&minusOne, sol, errorCtx->solution, e);                                                 CHKERRQ(ierr);
  ierr = GVecView(e, errorCtx->error_viewer);                                                             CHKERRQ(ierr);
  ierr = GTSDestroyContext(ts, (PetscObject) errorCtx->ctx, oldCtx);                                      CHKERRQ(ierr);

  /* Compute 2-norm and max-norm of error */
  if (errorCtx->norm_error_viewer) {
    ierr = PetscObjectGetComm((PetscObject) ts, &comm);                                                   CHKERRQ(ierr);
    ierr = PetscViewerASCIIGetPointer(errorCtx->norm_error_viewer, &file);                                CHKERRQ(ierr);
    ierr = VecNorm(e, NORM_2, &norm_2);                                                                   CHKERRQ(ierr);
    ierr = VecNorm(e, NORM_MAX, &norm_max);                                                               CHKERRQ(ierr);
    ierr = VecNorm(errorCtx->solution, NORM_2, &sol_norm_2);                                              CHKERRQ(ierr);
    PetscFPrintf(comm, file, "Timestep %d time %g error 2-norm %g error max-norm %g exact sol 2-norm %g\n",
                 step, ltime, norm_2, norm_max, sol_norm_2);
  }
  ierr = GVecRestoreWorkGVec(errorCtx->solution, &e);                                                     CHKERRQ(ierr);  
  PetscFunctionReturn(0);
}

EXTERN_C_BEGIN
#undef  __FUNCT__
#define __FUNCT__ "GTSOptionsChecker_Private"
int GTSOptionsChecker_Private(GTS ts) 
{
  char       *prefix;  
  int         loc[4], nmax;
  MPI_Comm    comm;
  PetscViewer viewer;
  PetscDraw   draw;
  PetscTruth  opt;
  int         ierr;

  PetscFunctionBegin;
  ierr = TSGetOptionsPrefix(ts, &prefix);                                                                 CHKERRQ(ierr);
  ierr = PetscObjectGetComm((PetscObject) ts, &comm);                                                     CHKERRQ(ierr);

  nmax = 4;
  loc[0] = PETSC_DECIDE; loc[1] = PETSC_DECIDE; loc[2] = 300; loc[3] = 300;
  ierr = PetscOptionsGetIntArray(prefix, "-gvec_ts_solutionmonitor", loc, &nmax, &opt);                   CHKERRQ(ierr);
  if (opt) {
    ierr = PetscViewerDrawOpen(comm, 0, 0, loc[0], loc[1], loc[2], loc[3], &viewer);                      CHKERRQ(ierr);
    ierr = PetscViewerDrawGetDraw(viewer, 0, &draw);                                                      CHKERRQ(ierr);
    ierr = PetscDrawSetTitle(draw, "Approx. Solution");                                                   CHKERRQ(ierr);
    PetscLogObjectParent(ts, (PetscObject) viewer);
    ierr = TSSetMonitor(ts, GTSSolutionMonitor, (void *) viewer, PETSC_NULL);                             CHKERRQ(ierr);
  }
#if 0
  nmax = 4;
  loc[0] = PETSC_DECIDE; loc[1] = PETSC_DECIDE; loc[2] = 300; loc[3] = 300;
  ierr = PetscOptionsGetIntArray(prefix, "-gvec_ts_errormonitor", loc, &nmax, &opt);                      CHKERRQ(ierr);
  if (opt) {
    ierr = PetscViewerDrawOpenX(comm, 0, 0, loc[0], loc[1], loc[2], loc[3], &viewer);                     CHKERRQ(ierr);
    ierr = PetscViewerDrawGetDraw(viewer, &draw);                                                         CHKERRQ(ierr);
    ierr = PetscDrawSetTitle(draw, "Error");                                                              CHKERRQ(ierr);
    PetscLogObjectParent(ts, (PetscObject) viewer);
    ierr = TSSetMonitor(ts, GTSErrorMonitor, (void *) viewer);                                            CHKERRQ(ierr);
  }
#endif
  /*--------------------------------------------------------------------- */
  ierr = PetscOptionsHasName(PETSC_NULL, "-help", &opt);                                                  CHKERRQ(ierr);
  if (opt) {
    char *pprefix;
    int   len = 2;
    if (prefix != PETSC_NULL) {
      ierr = PetscStrlen(prefix, &len);                                                                   CHKERRQ(ierr);
    }
    ierr = PetscMalloc((len+2) * sizeof(char), &pprefix);                                                 CHKERRQ(ierr);
    PetscStrcpy(pprefix, "-");
    if (prefix != PETSC_NULL)
      PetscStrcat(pprefix, prefix);
    PetscPrintf(comm," Additional TS Monitor options for grid vectors\n");
    PetscPrintf(comm,"   %sgvec_ts_solutionmonitor\n", pprefix);
    PetscPrintf(comm,"   %sgvec_ts_errormonitor -- not working yet\n", pprefix);
    ierr = PetscFree(pprefix);                                                                            CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}
EXTERN_C_END

/*---------------------------------------------- Context Functions --------------------------------------------------*/
EXTERN_C_BEGIN
#undef  __FUNCT__
#define __FUNCT__ "GTSContextDestroy_Private"
int GTSContextDestroy_Private(GTSContext ctx)
{
  PetscFunctionBegin;
  if (--ctx->refct > 0)
    PetscFunctionReturn(0);
  PetscLogObjectDestroy(ctx);
  PetscHeaderDestroy(ctx);
  PetscFunctionReturn(0);
}
EXTERN_C_END

#undef  __FUNCT__
#define __FUNCT__ "GTSCreateContext"
/*@
  GTSCreateContext - This functions creates a GTSContext attached to the given context.

  Collective on GTS

  Input Parameters:
. ts      - The GTS
. ltime   - The time for the context
. ctx     - The parent context

  Output Parameter:
. oldCtx  - The previous GTSContext

  Level: developer

.keywords: grid ts, context, create
.seealso: GTSDestroyContext(), GTSCreateConstraintContext()
@*/
int GTSCreateContext(GTS ts, double ltime, PetscObject ctx, PetscObject *oldCtx)
{
  GTSContext GTSctx;
  int        ierr;

  PetscFunctionBegin;
  /* Preserve previous context */
  ierr = PetscObjectQuery(ctx, "GTSContext", oldCtx);                                                     CHKERRQ(ierr);
  if (*oldCtx != PETSC_NULL) {
    ierr = PetscObjectReference(*oldCtx);                                                                 CHKERRQ(ierr);
  }
  /* Setup new context */
  PetscHeaderCreate(GTSctx, _GTSContext, int, TS_COOKIE, -1, "context", ts->comm, GTSContextDestroy_Private, 0);
  PetscLogObjectCreate(GTSctx);
  GTSctx->t     = ltime;
  GTSctx->dt    = ts->time_step;
  GTSctx->sol   = ts->vec_sol;
  GTSctx->nwork = ts->nwork;
  GTSctx->work  = ts->work;
  ierr = PetscObjectCompose(ctx, "GTSContext", (PetscObject) GTSctx);                                     CHKERRQ(ierr);
  /* This makes sure the context is destroyed when the object is destroyed */
  ierr = PetscObjectDereference((PetscObject) GTSctx);                                                    CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "GTSDestroyContext"
/*@
  GTSDestroyContext - This functions destroys a GTSContext attached to the given context and
  reattaches the old context.

  Collective on GTS

  Input Parameters:
. ts     - The GTS
. ctx    - The parent context
. oldCtx - The previous GTSContext

  Level: developer

.keywords: grid ts, context, destroy
.seealso: GTSCreateContext(), GTSCreateConstraintContext()
@*/
int GTSDestroyContext(GTS ts, PetscObject ctx, PetscObject oldCtx)
{
  int ierr;

  PetscFunctionBegin;
  if (oldCtx != PETSC_NULL) {
    /* Restore previous context -- New one is destroyed automatically */
    ierr = PetscObjectCompose(ctx, "GTSContext", oldCtx);                                                 CHKERRQ(ierr);
    ierr = PetscObjectDereference(oldCtx);                                                                CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "GTSCreateConstraintContext"
/*@
  GTSCreateConstraintContext - This functions creates a GTSContext attached to the
  constraint context of the underlying grid.

  Collective on GTS

  Input Parameter:
. ts  - The GTS

  Level: developer

.keywords: grid ts, context
.seealso: GTSCreateContext(), GTSDestroyContext()
@*/
int GTSCreateConstraintContext(GTS ts)
{
  Grid                  grid;
  PetscTruth            isConstrained;
  PetscConstraintObject constCtx;
  PetscObject           oldCtx;
  int                   ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ts, TS_COOKIE);
  if (ts->isGTS == PETSC_FALSE)
    PetscFunctionReturn(0);
  ierr = GTSGetGrid(ts, &grid);                                                                           CHKERRQ(ierr);
  ierr = GridIsConstrained(grid, &isConstrained);                                                         CHKERRQ(ierr);
  if (isConstrained == PETSC_TRUE) {
    ierr = GridGetConstraintContext(grid, &constCtx);                                                     CHKERRQ(ierr);
    ierr = GTSCreateContext(ts, ts->ptime, (PetscObject) constCtx, &oldCtx);                              CHKERRQ(ierr);
    /* Remove the extra reference from GTSCreateContext() */
    if (oldCtx != PETSC_NULL) {
      ierr = PetscObjectDereference(oldCtx);                                                              CHKERRQ(ierr);
    }
  }
  PetscFunctionReturn(0);
}