#ifdef PETSC_RCS_HEADER
static char vcid[] = "$Id: gbeuler.c,v 1.24 2000/01/10 03:18:55 knepley Exp $";
#endif
/*
       Code for Timestepping with implicit backwards Euler using Grid Solvers.
*/
#include "src/ts/tsimpl.h"                /*I   "ts.h"   I*/
#include "gsolver.h"

typedef struct {
  Vec update;      /* work vector where new solution is formed */
  Vec func;        /* work vector where F(t[i],u[i]) is stored */
  Vec rhs;         /* work vector for RHS; vec_sol/dt */
} TS_BEuler;

extern int TSCreate_BEuler(TS);

#undef __FUNCT__  
#define __FUNCT__ "GTSStep_BEuler_Linear"
static int GTSStep_BEuler_Linear(GTS ts, int *steps, PetscReal *ltime)
{
  TS_BEuler    *beuler    = (TS_BEuler*) ts->data;
  int           max_steps = ts->max_steps;
  Grid          grid;
  MatStructure  str;
  PetscTruth    isTimeDependent;
  PetscScalar   mdt;
  int           numFields, field, its;
  int           step, f;
  int           ierr;

  PetscFunctionBegin;
  *steps = -ts->steps;
  ierr   = TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol);                                              CHKERRQ(ierr);

  /* Set initial guess to be previous solution */
  ierr = VecCopy(ts->vec_sol, beuler->update);                                                            CHKERRQ(ierr);

  /* Step through time */
  ierr = GTSGetGrid(ts, &grid);                                                                           CHKERRQ(ierr);
  ierr = GridGetNumActiveFields(grid, &numFields);                                                        CHKERRQ(ierr);
  for(step = 0; step < max_steps; step++) {
    /* Update mesh */
    ts->time_step_old = ts->time_step;
    ierr = (*ts->ops->update)(ts, ts->ptime, &ts->time_step);                                             CHKERRQ(ierr);
    if (ts->time_step != ts->time_step_old) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Variable time step not allowed");

    /* Increment time */
    ts->ptime += ts->time_step;
    mdt        = 1.0/ts->time_step;
    if (ts->ptime > ts->max_time) break;

    /* Update the boundary values */
    ierr = GTSCalcBCValues(ts);                                                                           CHKERRQ(ierr);

    /* Compute A(t) and B(t) */
    ierr = GTSEvaluateSystemMatrix(ts, ts->ptime, &ts->A, &ts->B, &str, (PetscObject) ts->jacP);          CHKERRQ(ierr);

    /* Compute f + \Delta t^{-1} I U^n */
    ierr = GTSEvaluateRhs(ts, ts->ptime, ts->vec_sol, beuler->rhs, (PetscObject) ts->funP);               CHKERRQ(ierr);
    for(f = 0; f < numFields; f++) {
      ierr = GridGetActiveField(grid, f, &field);                                                         CHKERRQ(ierr);
      ierr = GTSGetTimeDependence(ts, field, &isTimeDependent);                                           CHKERRQ(ierr);
      if (isTimeDependent == PETSC_TRUE) {
        ierr = GVecEvaluateOperatorGalerkin(beuler->rhs, ts->vec_sol, ts->vec_sol, field, field, IDENTITY, mdt, ts->funP); CHKERRQ(ierr);
      }
    }
    ierr = (*ts->ops->applyrhsbc)(ts, beuler->rhs, ts->funP);                                             CHKERRQ(ierr);

    /* Solve (\Delta t^{-1} I + A) U^{n+1} = \Delta t^{-1} I U^n */
    ierr = SLESSetOperators(ts->sles, ts->A, ts->B, str);                                                 CHKERRQ(ierr);
    ierr = SLESSolve(ts->sles, beuler->rhs, beuler->update, &its);                                        CHKERRQ(ierr);
    ts->linear_its += PetscAbsInt(its);
    ts->steps++;

    /* Update solution U^n --> U^{n+1} */
    ierr = VecCopy(beuler->update, ts->vec_sol);                                                          CHKERRQ(ierr);
    ierr = TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol);                                              CHKERRQ(ierr);
  }

  *steps  += ts->steps;
  *ltime   = ts->ptime;
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "GTSStep_BEuler_Nonlinear"
static int GTSStep_BEuler_Nonlinear(GTS ts, int *steps, PetscReal *ltime)
{
  TS_BEuler *beuler    = (TS_BEuler*) ts->data;
  Grid       grid;
  int        max_steps = ts->max_steps;
  int        its, lits, numFields;
  int        step, field, f;
  int        ierr;
  
  PetscFunctionBegin;
  *steps = -ts->steps;
  ierr   = TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol);                                             CHKERRQ(ierr);

  /* Set initial guess to be previous solution */
  ierr = VecCopy(ts->vec_sol, beuler->update);                                                           CHKERRQ(ierr);
  ierr = GTSGetGrid(ts, &grid);                                                                          CHKERRQ(ierr);
  ierr = GridGetNumActiveFields(grid, &numFields);                                                       CHKERRQ(ierr);

  /* Step through time */
  for(step = 0; step < max_steps; step++)
  {
    /* Update mesh */
    ts->time_step_old = ts->time_step;
    ierr = (*ts->ops->update)(ts, ts->ptime, &ts->time_step);                                            CHKERRQ(ierr);
    if (ts->time_step != ts->time_step_old)
    {
      /* Multiply \Delta t^{-1}_0 I by {\Delta t_0 \over \Delta t} */
      for(f = 0; f < numFields; f++)
      {
        ierr = GridGetActiveField(grid, f, &field);                                                      CHKERRQ(ierr);
        if (ts->Iindex[field] >= 0)
          {ierr = GridScaleMatOperator(grid, ts->time_step_old/ts->time_step, ts->Iindex[field]);        CHKERRQ(ierr);}
      }
      ts->time_step_old = ts->time_step;
    }

    /* Increment time */
    ts->ptime += ts->time_step;
    if (ts->ptime > ts->max_time) break;

    /* Update the boundary values */
    ierr = GTSCalcBCValues(ts);                                                                          CHKERRQ(ierr);

    /* Solve \Delta t^{-1} I (U^{n+1} - U^n) - F(U^{n+1}, t^{n+1}) = 0 */
    ierr = SNESSolve(ts->snes, beuler->update, &its);                                                    CHKERRQ(ierr);
    ierr = SNESGetNumberLinearIterations(ts->snes, &lits);                                               CHKERRQ(ierr);
    ts->nonlinear_its += PetscAbsInt(its);
    ts->linear_its    += PetscAbsInt(lits);
    ts->steps++;

    /* Update solution U^n --> U^{n+1} */
    ierr = VecCopy(beuler->update, ts->vec_sol);                                                         CHKERRQ(ierr);
    ierr = TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol);                                             CHKERRQ(ierr);
  }

  *steps  += ts->steps;
  *ltime   = ts->ptime;
  PetscFunctionReturn(0);
}

/* 
   This defines the nonlinear equation that is to be solved with GSNES

     \Delta t^{-1} I (U^{n+1} - U^{n}) - F(U^{n+1}, t^{n+1})
*/
#undef __FUNCT__  
#define __FUNCT__ "GTSBEulerFunction"
int GTSBEulerFunction(GSNES snes, GVec x, GVec y, void *ctx)
{
  Grid        grid;
  PetscObject oldCtx;
  TS          ts;
  PetscTruth  isTimeDependent;
  PetscScalar mone = -1.0;
  PetscScalar mdt;
  int         numFields, field, f;
  int         ierr;

  PetscFunctionBegin;
  ierr = PetscObjectQuery((PetscObject) ctx, "TS", (PetscObject *) &ts);                                  CHKERRQ(ierr);
  PetscValidHeaderSpecific(ts, TS_COOKIE);
  mdt  = 1.0/ts->time_step;
  /* Make -F(U^{n+1}, t^{n+1}) */
  ierr = GTSEvaluateRhs(ts, ts->ptime, x, y, (PetscObject) ts->funP);                                     CHKERRQ(ierr);

  /* Add \Delta t^{-1} I (U^{n+1} - U^n) for each time dependent field */
  ierr = GTSGetGrid(ts, &grid);                                                                           CHKERRQ(ierr);
  ierr = GTSCreateContext(ts, ts->ptime, (PetscObject) ts->funP, &oldCtx);                                CHKERRQ(ierr);
  ierr = VecWAXPY(&mone, ts->vec_sol, x, ts->work[0]);                                                    CHKERRQ(ierr);
  ierr = GridCalcBCValuesDifference(grid);                                                                CHKERRQ(ierr);
  ierr = GridSetBCValuesType(grid, BC_VALUES_DIFF);                                                       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 = GVecEvaluateOperatorGalerkin(y, ts->work[0], ts->work[0], field, field, IDENTITY, mdt, ts->funP); CHKERRQ(ierr);
    }
  }
  ierr = GridSetBCValuesType(grid, BC_VALUES);                                                            CHKERRQ(ierr);

  /* Apply boundary conditions */
  ierr = (*ts->ops->applyrhsbc)(ts, y, ts->funP);                                                         CHKERRQ(ierr);

  ierr = GTSDestroyContext(ts, (PetscObject) ts->funP, oldCtx);                                           CHKERRQ(ierr);
  PetscFunctionReturn(0);
}
/* 
   This defines the nonlinear equation that is to be solved with GSNES

     P^T (\Delta t^{-1} I P (U^{n+1} - U^{n}) - F(P U^{n+1}, t^{n+1})) + \Delta t^{-1} I_P (U^{n+1}_P - U^n_P) - F_g
*/
#undef __FUNCT__  
#define __FUNCT__ "GTSBEulerFunctionConstrained"
int GTSBEulerFunctionConstrained(GSNES snes, GVec x, GVec y, void *ctx)
{
  Grid        grid;
  PetscObject oldCtx;
  TS          ts;
  Mat         P;
  PetscTruth  isTimeDependent;
  PetscScalar mdt;
  PetscScalar mone = -1.0;
  int         numFields, field, f;
  int         ierr;

  PetscFunctionBegin;
  ierr = PetscObjectQuery((PetscObject) ctx, "TS", (PetscObject *) &ts);                                  CHKERRQ(ierr);
  PetscValidHeaderSpecific(ts, TS_COOKIE);
  mdt  = 1.0/ts->time_step;
#ifdef PETSC_USE_BOPT_g
  ierr = PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);                                           CHKERRQ(ierr);
#endif

  /* Make -F(P U^{n+1}, t^{n+1}) */
  ierr = GTSGetGrid(ts, &grid);                                                                           CHKERRQ(ierr);
  ierr = GTSCreateContext(ts, ts->ptime, (PetscObject) ts->funP, &oldCtx);                                CHKERRQ(ierr);
  ierr = GridGetConstraintMatrix(grid, &P);                                                               CHKERRQ(ierr);
  ierr = MatMult(P, x, ts->work[0]);                                                                      CHKERRQ(ierr);
  ierr = GTSEvaluateRhs(ts, ts->ptime, ts->work[0], ts->work[1], (PetscObject) ts->funP);                 CHKERRQ(ierr);

  /* Add \Delta t^{-1} I P (U^{n+1} - U^n) for each time dependent field */
  ierr = MatMult(P, ts->vec_sol, ts->work[2]);                                                            CHKERRQ(ierr);
  ierr = VecWAXPY(&mone, ts->work[2], ts->work[0], ts->work[0]);                                          CHKERRQ(ierr);
  ierr = GridCalcBCValuesDifference(grid);                                                                CHKERRQ(ierr);
  ierr = GridSetBCValuesType(grid, BC_VALUES_DIFF);                                                       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 = GVecEvaluateOperatorGalerkin(ts->work[1], ts->work[0], ts->work[0], field, field, IDENTITY, mdt, ts->funP); CHKERRQ(ierr);
    }
  }
  ierr = GridSetBCValuesType(grid, BC_VALUES);                                                            CHKERRQ(ierr);

  /* Reduce the system with constraints: apply P^T  */
  ierr = MatMultTranspose(P, ts->work[1], y);                                                             CHKERRQ(ierr);

  /* Add in extra degrees of freedom */
  ierr = (*((PetscConstraintObject) ctx)->ops->applyrhs)(grid, x, y);                                     CHKERRQ(ierr);

  /* Apply boundary conditions */
  ierr = (*ts->ops->applyrhsbc)(ts, y, ts->funP);                                                         CHKERRQ(ierr);

  ierr = GTSDestroyContext(ts, (PetscObject) ts->funP, oldCtx);                                           CHKERRQ(ierr);
#ifdef PETSC_USE_BOPT_g
  ierr = PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);                                           CHKERRQ(ierr);
#endif
  PetscFunctionReturn(0);
}

/*
   This constructs the Jacobian needed for GSNES 

     J = \Delta t^{-1} I - J_{F}

   where J_{F} is the given Jacobian of F.
*/
#undef __FUNCT__  
#define __FUNCT__ "GTSBEulerJacobian"
int GTSBEulerJacobian(GSNES snes, GVec x, GMat *J, GMat *M, MatStructure *str, void *ctx)
{
  GTS ts;
  int ierr;

  PetscFunctionBegin;
  ierr = PetscObjectQuery((PetscObject) ctx, "TS", (PetscObject *) &ts);                                  CHKERRQ(ierr);
  PetscValidHeaderSpecific(ts, TS_COOKIE);
  ierr = GTSEvaluateJacobian(ts, ts->ptime, x, J, M, str, (PetscObject) ctx);                             CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*
   This constructs the Jacobian needed for GSNES 

     J = \Delta t^{-1} I - J_{F}(P U^n)

   where J_{F} is the given Jacobian of F.
*/
#undef __FUNCT__  
#define __FUNCT__ "GTSBEulerJacobianConstrained"
int GTSBEulerJacobianConstrained(GSNES snes, GVec x, GMat *J, GMat *M, MatStructure *str, void *ctx)
{
  GTS  ts;
  Grid grid;
  Mat  P;
  int  ierr;

  PetscFunctionBegin;
  ierr = PetscObjectQuery((PetscObject) ctx, "TS", (PetscObject *) &ts);                                  CHKERRQ(ierr);
  PetscValidHeaderSpecific(ts, TS_COOKIE);
  ierr = GTSGetGrid(ts, &grid);                                                                           CHKERRQ(ierr);
  ierr = GridGetConstraintMatrix(grid, &P);                                                               CHKERRQ(ierr);
  ierr = MatMult(P, x, ts->work[0]);                                                                      CHKERRQ(ierr);
  ierr = GTSEvaluateJacobian(ts, ts->ptime, ts->work[0], J, M, str, (PetscObject) ctx);                   CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "GTSSetUp_BEuler_Linear"
static int GTSSetUp_BEuler_Linear(GTS ts) {
  TS_BEuler *beuler = (TS_BEuler*) ts->data;
  Grid       grid;
  int        ierr;

  PetscFunctionBegin;
  ierr = GTSGetGrid(ts, &grid);                                                                          CHKERRQ(ierr);
  ierr = VecDuplicate(ts->vec_sol, &beuler->update);                                                     CHKERRQ(ierr);
  ierr = VecDuplicate(ts->vec_sol, &beuler->rhs);                                                        CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "GTSSetupGSNES_BEuler"
static int GTSSetupGSNES_BEuler(GTS ts)
{
  GVec func;
  GMat A, B;
  int  ierr;

  PetscFunctionBegin;
  ierr = SNESGetFunction(ts->snes, &func, PETSC_NULL, PETSC_NULL);                                        CHKERRQ(ierr);
  ierr = PetscObjectCompose((PetscObject) ts->funP, "TS", (PetscObject) ts);                              CHKERRQ(ierr);
  ierr = SNESSetFunction(ts->snes, func, GTSBEulerFunction, ts->funP);                                    CHKERRQ(ierr);
  ierr = SNESGetJacobian(ts->snes, &A, &B, PETSC_NULL, PETSC_NULL);                                       CHKERRQ(ierr);
  ierr = PetscObjectCompose((PetscObject) ts->jacP, "TS", (PetscObject) ts);                              CHKERRQ(ierr);
  ierr = SNESSetJacobian(ts->snes, A, B, GTSBEulerJacobian, ts->jacP);                                    CHKERRQ(ierr);

  /* This requires that ts is the context for the function */
  ierr = SNESSetSolutionBC(ts->snes, GTSSolutionBCforGSNES);                                              CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "GTSSetUp_BEuler_Nonlinear"
static int GTSSetUp_BEuler_Nonlinear(GTS ts)
{
  TS_BEuler *beuler = (TS_BEuler*) ts->data;
  Grid       grid;
  int        ierr;

  PetscFunctionBegin;
  ierr = GTSGetGrid(ts, &grid);                                                                          CHKERRQ(ierr);

  /* Setup Rhs \Delta t^{-1} I (U^{n+1} - U^{n}) - F(U^{n+1}, t^{n+1}) */
  ierr = VecDuplicate(ts->vec_sol, &beuler->update);                                                     CHKERRQ(ierr);
  ts->nwork = 1;
  ierr = VecDuplicateVecs(ts->vec_sol, ts->nwork, &ts->work);                                            CHKERRQ(ierr);

  /* Setup the nonlinear solver */
  ierr = GTSSetupGSNES_BEuler(ts);                                                                       CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "GTSSetupGSNES_BEuler"
static int GTSSetupGSNES_BEuler_Constrained(GTS ts)
{
  GVec func;
  GMat A, B;
  int  ierr;

  PetscFunctionBegin;
  ierr = SNESGetFunction(ts->snes, &func, PETSC_NULL, PETSC_NULL);                                        CHKERRQ(ierr);
  ierr = PetscObjectCompose((PetscObject) ts->funP, "TS", (PetscObject) ts);                              CHKERRQ(ierr);
  ierr = SNESSetFunction(ts->snes, func, GTSBEulerFunctionConstrained, ts->funP);                         CHKERRQ(ierr);
  ierr = SNESGetJacobian(ts->snes, &A, &B, PETSC_NULL, PETSC_NULL);                                       CHKERRQ(ierr);
  ierr = PetscObjectCompose((PetscObject) ts->jacP, "TS", (PetscObject) ts);                              CHKERRQ(ierr);
  ierr = SNESSetJacobian(ts->snes, A, B, GTSBEulerJacobianConstrained, ts->jacP);                         CHKERRQ(ierr);

  /* This requires that ts is the context for the function */
  ierr = SNESSetSolutionBC(ts->snes, GTSSolutionBCforGSNES);                                              CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "GTSSetUp_BEuler_Nonlinear_Constrained"
static int GTSSetUp_BEuler_Nonlinear_Constrained(GTS ts)
{
  TS_BEuler  *beuler = (TS_BEuler*) ts->data;
  Grid        grid;
  PetscTruth  explicitConst;
  int         ierr;

  PetscFunctionBegin;
  ierr = GTSGetGrid(ts, &grid);                                                                           CHKERRQ(ierr);

  /* Setup Rhs \Delta t^{-1} I (U^{n+1} - U^{n}) - F(U^{n+1}, t^{n+1}) */
  ierr = VecDuplicate(ts->vec_sol, &beuler->update);                                                      CHKERRQ(ierr);
  ts->nwork = 3;
  ierr = PetscMalloc(ts->nwork * sizeof(GVec *), &ts->work);                                              CHKERRQ(ierr);

  /* Setup the nonlinear solver */
  ierr = GridGetExplicitConstraints(grid, &explicitConst);                                                CHKERRQ(ierr);
  if (explicitConst == PETSC_FALSE) {
    ierr = GTSSetupGSNES_BEuler_Constrained(ts);                                                          CHKERRQ(ierr);
    ierr = GVecCreate(grid, &ts->work[0]);                                                                CHKERRQ(ierr);
  } else {
    ierr = GTSSetupGSNES_BEuler(ts);                                                                      CHKERRQ(ierr);
    ierr = GVecCreateConstrained(grid, &ts->work[0]);                                                     CHKERRQ(ierr);
  }
  ierr = GVecDuplicate(ts->work[0], &ts->work[1]);                                                        CHKERRQ(ierr);
  ierr = GVecDuplicate(ts->work[0], &ts->work[2]);                                                        CHKERRQ(ierr);

  /* Set the context */
  ierr = GTSCreateConstraintContext(ts);                                                                  CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "GTSReform_BEuler"
static int GTSReform_BEuler(GTS ts)
{
  PetscFunctionBegin;
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "GTSReallocate_BEuler"
static int GTSReallocate_BEuler(GTS ts)
{
  TS_BEuler *beuler = (TS_BEuler*) ts->data;
  Grid       grid;
  GSNES      newSnes;
  int        ierr;

  PetscFunctionBegin;
  /* Destroy old structures */
  ierr = VecDestroy(beuler->update);                                                                     CHKERRQ(ierr);
  if (ts->nwork) {
    ierr = VecDestroyVecs(ts->work, ts->nwork);                                                          CHKERRQ(ierr);
  }

  /* Recreate GSNES */
  if (ts->snes) {
    ierr = GTSGetGrid(ts, &grid);                                                                        CHKERRQ(ierr);
    ierr = GSNESCreate(grid, ts, &newSnes);                                                              CHKERRQ(ierr);
    ierr = SNESSetFromOptions(newSnes);                                                                  CHKERRQ(ierr);
    ierr = GSNESDuplicateMonitors(ts->snes, newSnes);                                                    CHKERRQ(ierr);
    ierr = GSNESDestroy(ts->snes);                                                                       CHKERRQ(ierr);
    ts->snes = newSnes;
  }

  /* Recreate structures */
  ierr = TSSetUp(ts);                                                                                    CHKERRQ(ierr);
    
  /* Start off with the current solution */
  ierr = VecCopy(ts->vec_sol, beuler->update);                                                           CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "GTSDestroy_BEuler"
static int GTSDestroy_BEuler(TS ts)
{
  TS_BEuler *beuler = (TS_BEuler*) ts->data;
  int        ierr;

  PetscFunctionBegin;
  if (beuler->update) {
    ierr = VecDestroy(beuler->update);                                                                    CHKERRQ(ierr);
  }
  if (beuler->func) {
    ierr = VecDestroy(beuler->func);                                                                      CHKERRQ(ierr);
  }
  if (beuler->rhs) {
    ierr = VecDestroy(beuler->rhs);                                                                       CHKERRQ(ierr);
  }
  if (ts->nwork) {
    ierr = VecDestroyVecs(ts->work, ts->nwork);                                                           CHKERRQ(ierr);
  }
  ierr = PetscFree(beuler);                                                                               CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "GTSSetFromOptions_BEuler_Linear"
static int GTSSetFromOptions_BEuler_Linear(TS ts)
{
  PetscFunctionBegin;
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "GTSSetFromOptions_BEuler_Nonlinear"
static int GTSSetFromOptions_BEuler_Nonlinear(TS ts)
{
  PetscFunctionBegin;
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "GTSPrintHelp_BEuler"
static int GTSPrintHelp_BEuler(TS ts, char *prefix)
{
  PetscFunctionBegin;
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "GTSView_BEuler"
static int GTSView_BEuler(TS ts, PetscViewer viewer)
{
  PetscFunctionBegin;
  PetscFunctionReturn(0);
}

EXTERN_C_BEGIN
#undef __FUNCT__  
#define __FUNCT__ "GTSCreate_BEuler"
int GTSCreate_BEuler(TS ts)
{
  TS_BEuler *beuler;
  KSP        ksp;
  Grid       grid;
  PetscTruth flag;
  int        ierr;

  PetscFunctionBegin;
  ts->ops->destroy   = GTSDestroy_BEuler;
  ts->ops->view      = GTSView_BEuler;
  ts->ops->printhelp = GTSPrintHelp_BEuler;
  ierr = PetscObjectChangeSerializeName((PetscObject) ts, GTS_SER_BEULER_BINARY);                         CHKERRQ(ierr);

  if (ts->problem_type == TS_LINEAR) {
    if (!ts->A) SETERRQ(PETSC_ERR_ARG_WRONG, "Must set rhs matrix for linear problem");
    ts->ops->setup          = GTSSetUp_BEuler_Linear;
    ts->ops->step           = GTSStep_BEuler_Linear;
    ts->ops->setfromoptions = GTSSetFromOptions_BEuler_Linear;
    ts->ops->reallocate     = GTSReallocate_BEuler;
    ts->ops->reform         = GTSReform_BEuler;
    ierr = SLESCreate(ts->comm, &ts->sles);                                                               CHKERRQ(ierr);
    ierr = SLESGetKSP(ts->sles, &ksp);                                                                    CHKERRQ(ierr);
    ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);                                                    CHKERRQ(ierr);
    ierr = GTSGetGrid(ts, &grid);                                                                         CHKERRQ(ierr);
    ierr = GridIsConstrained(grid, &flag);                                                                CHKERRQ(ierr);
    if (flag == PETSC_TRUE) SETERRQ(PETSC_ERR_SUP, "Constraints not supported for the linear problem yet");
  } else if (ts->problem_type == TS_NONLINEAR) {
    ierr = GTSGetGrid(ts, &grid);                                                                         CHKERRQ(ierr);
    ierr = GridIsConstrained(grid, &flag);                                                                CHKERRQ(ierr);
    ierr = GSNESCreate(grid, ts, &ts->snes);                                                              CHKERRQ(ierr);
    if (flag) {
      ts->ops->setup          = GTSSetUp_BEuler_Nonlinear_Constrained;  
      ts->ops->step           = GTSStep_BEuler_Nonlinear;
      ts->ops->setfromoptions = GTSSetFromOptions_BEuler_Nonlinear;
      ts->ops->reform         = GTSReform_BEuler;
      ts->ops->reallocate     = GTSReallocate_BEuler;
    } else {
      ts->ops->setup          = GTSSetUp_BEuler_Nonlinear;  
      ts->ops->step           = GTSStep_BEuler_Nonlinear;
      ts->ops->setfromoptions = GTSSetFromOptions_BEuler_Nonlinear;
      ts->ops->reform         = GTSReform_BEuler;
      ts->ops->reallocate     = GTSReallocate_BEuler;
    }
  } else {
    SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid problem type");
  }

  ierr = PetscNew(TS_BEuler, &beuler);                                                                    CHKERRQ(ierr);
  PetscLogObjectMemory(ts, sizeof(TS_BEuler));
  ierr = PetscMemzero(beuler, sizeof(TS_BEuler));                                                         CHKERRQ(ierr);
  ts->data = (void *) beuler;

  PetscFunctionReturn(0);
}
EXTERN_C_END

EXTERN_C_BEGIN
#undef __FUNCT__  
#define __FUNCT__ "GTSSerialize_BEuler"
int GTSSerialize_BEuler(MPI_Comm comm, TS *ts, PetscViewer viewer, PetscTruth store)
{
  TS   t;
  Grid grid;
  int  fd;
  int  numFields;
  int  hasPC;
  int  zero = 0;
  int  one  = 0;
  int  ierr;

  PetscFunctionBegin;
  ierr = PetscViewerBinaryGetDescriptor(viewer, &fd);                                                     CHKERRQ(ierr);
  if (store) {
    t = *ts;
    ierr = PetscBinaryWrite(fd, &t->problem_type,      1,         PETSC_INT,     0);                      CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd, &t->isGTS,             1,         PETSC_INT,     0);                      CHKERRQ(ierr);
    ierr = GTSGetGrid(t, &grid);                                                                          CHKERRQ(ierr);
    ierr = GridGetNumFields(grid, &numFields);                                                            CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd, t->isExplicit,         numFields, PETSC_INT,     0);                      CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd, t->Iindex,             numFields, PETSC_INT,     0);                      CHKERRQ(ierr);
    if (t->problem_type == TS_LINEAR) {
      ierr = GMatSerialize(grid, &t->A,    viewer, store);                                                CHKERRQ(ierr);
      if (t->B == t->A) {
        ierr = PetscBinaryWrite(fd, &zero,           1,         PETSC_INT,     0);                        CHKERRQ(ierr);
      } else {
        ierr = PetscBinaryWrite(fd, &one,            1,         PETSC_INT,     0);                        CHKERRQ(ierr);
        ierr = GMatSerialize(grid, &t->B,  viewer, store);                                                CHKERRQ(ierr);
      }
     }
    ierr = PetscBinaryWrite(fd, &t->max_steps,         1,         PETSC_INT,     0);                      CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd, &t->max_time,          1,         PETSC_SCALAR,  0);                      CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd, &t->time_step,         1,         PETSC_SCALAR,  0);                      CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd, &t->initial_time_step, 1,         PETSC_SCALAR,  0);                      CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd, &t->max_time,          1,         PETSC_SCALAR,  0);                      CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd, &t->steps,             1,         PETSC_INT,     0);                      CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd, &t->ptime,             1,         PETSC_SCALAR,  0);                      CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd, &t->linear_its,        1,         PETSC_INT,     0);                      CHKERRQ(ierr);
    ierr = PetscBinaryWrite(fd, &t->nonlinear_its,     1,         PETSC_INT,     0);                      CHKERRQ(ierr);

    /* Protect against nasty option overwrites */
    ierr = PetscOptionsClearValue("-ts_init_time");                                                       CHKERRQ(ierr);
    ierr = PetscOptionsClearValue("-ts_init_time_step");                                                  CHKERRQ(ierr);
  } else {
    ierr = TSCreate(comm, &t);                                                                            CHKERRQ(ierr);

    ierr = PetscBinaryRead(fd, &t->problem_type,      1,         PETSC_INT);                              CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd, &t->isGTS,             1,         PETSC_INT);                              CHKERRQ(ierr);
    t->bops->destroy = (int (*)(PetscObject)) GTSDestroy;
    t->bops->view    = (int (*)(PetscObject, PetscViewer)) GTSView;
    ierr = PetscStrallocpy("gbeuler", &t->type_name);                                                     CHKERRQ(ierr);

    grid = (Grid) *ts;
    ierr = GridGetNumFields(grid, &numFields);                                                            CHKERRQ(ierr);
    ierr = PetscMalloc(numFields * sizeof(PetscTruth), &t->isExplicit);                                   CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd, t->isExplicit,         numFields, PETSC_INT);                              CHKERRQ(ierr);
    ierr = PetscMalloc(numFields * sizeof(int),        &t->Iindex);                                       CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd, t->Iindex,             numFields, PETSC_INT);                              CHKERRQ(ierr);
    if (t->problem_type == TS_LINEAR) {
      ierr = GMatSerialize(grid, &t->A,    viewer, store);                                                CHKERRQ(ierr);
      ierr = PetscBinaryRead(fd, &hasPC,            1,         PETSC_INT);                                CHKERRQ(ierr);
      if (hasPC) {
        ierr = GMatSerialize(grid, &t->B, viewer, store);                                                 CHKERRQ(ierr);
      } else {
        t->B = t->A;
      }
    }
    /* We must have a Grid parent for the constructor */
    ierr = PetscObjectCompose((PetscObject) t, "Grid", (PetscObject) grid);                               CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd, &t->max_steps,         1,         PETSC_INT);                              CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd, &t->max_time,          1,         PETSC_SCALAR);                           CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd, &t->time_step,         1,         PETSC_SCALAR);                           CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd, &t->initial_time_step, 1,         PETSC_SCALAR);                           CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd, &t->max_time,          1,         PETSC_SCALAR);                           CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd, &t->steps,             1,         PETSC_INT);                              CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd, &t->ptime,             1,         PETSC_SCALAR);                           CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd, &t->linear_its,        1,         PETSC_INT);                              CHKERRQ(ierr);
    ierr = PetscBinaryRead(fd, &t->nonlinear_its,     1,         PETSC_INT);                              CHKERRQ(ierr);

    ierr = GTSCreate_BEuler(t);                                                                            CHKERRQ(ierr);
    *ts  = t;
  }

  PetscFunctionReturn(0);
}
EXTERN_C_END