static char help[] = "Solves heat equation in 1d.\n";

/*
  Solves the equation

    u_t = kappa  \Delta u
       Periodic boundary conditions

Evolve the  heat equation:
---------------
./heat -ts_monitor -snes_monitor  -pc_type lu  -draw_pause .1 -snes_converged_reason  -ts_type cn  -da_refine 5 -mymonitor

Evolve the  Allen-Cahn equation:
---------------
./heat -ts_monitor -snes_monitor  -pc_type lu  -draw_pause .1 -snes_converged_reason  -ts_type cn  -da_refine 5   -allen-cahn -kappa .001 -ts_max_time 5  -mymonitor

Evolve the  Allen-Cahn equation: zoom in on part of the domain
---------------
./heat -ts_monitor -snes_monitor  -pc_type lu   -snes_converged_reason     -ts_type cn  -da_refine 5   -allen-cahn -kappa .001 -ts_max_time 5  -zoom .25,.45  -mymonitor

The option -square_initial indicates it should use a square wave initial condition otherwise it loads the file InitialSolution.heat as the initial solution. You should run with
./heat -square_initial -ts_monitor -snes_monitor  -pc_type lu   -snes_converged_reason    -ts_type cn  -da_refine 9 -ts_max_time 1.e-4 -ts_time_step .125e-6 -snes_atol 1.e-25 -snes_rtol 1.e-25  -ts_max_steps 15
to generate InitialSolution.heat

*/
#include <petscdm.h>
#include <petscdmda.h>
#include <petscts.h>
#include <petscdraw.h>

/*
   User-defined routines
*/
extern PetscErrorCode FormFunction(TS, PetscReal, Vec, Vec, void *), FormInitialSolution(DM, Vec), MyMonitor(TS, PetscInt, PetscReal, Vec, void *), MyDestroy(void **);
typedef struct {
  PetscReal           kappa;
  PetscBool           allencahn;
  PetscDrawViewPorts *ports;
} UserCtx;

int main(int argc, char **argv)
{
  TS          ts;   /* time integrator */
  Vec         x, r; /* solution, residual vectors */
  PetscInt    steps, Mx;
  DM          da;
  PetscReal   dt;
  UserCtx     ctx;
  PetscBool   mymonitor;
  PetscViewer viewer;
  PetscBool   flg;

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Initialize program
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  PetscFunctionBeginUser;
  PetscCall(PetscInitialize(&argc, &argv, NULL, help));
  ctx.kappa = 1.0;
  PetscCall(PetscOptionsGetReal(NULL, NULL, "-kappa", &ctx.kappa, NULL));
  ctx.allencahn = PETSC_FALSE;
  PetscCall(PetscOptionsHasName(NULL, NULL, "-allen-cahn", &ctx.allencahn));
  PetscCall(PetscOptionsHasName(NULL, NULL, "-mymonitor", &mymonitor));

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Create distributed array (DMDA) to manage parallel grid and vectors
  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, 10, 1, 2, NULL, &da));
  PetscCall(DMSetFromOptions(da));
  PetscCall(DMSetUp(da));
  PetscCall(DMDASetFieldName(da, 0, "Heat equation: u"));
  PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
  dt = 1.0 / (ctx.kappa * Mx * Mx);

  /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Extract global vectors from DMDA; then duplicate for remaining
     vectors that are the same types
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  PetscCall(DMCreateGlobalVector(da, &x));
  PetscCall(VecDuplicate(x, &r));

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Create timestepping solver context
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
  PetscCall(TSSetDM(ts, da));
  PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
  PetscCall(TSSetRHSFunction(ts, NULL, FormFunction, &ctx));

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Customize nonlinear solver
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  PetscCall(TSSetType(ts, TSCN));

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Set initial conditions
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  PetscCall(FormInitialSolution(da, x));
  PetscCall(TSSetTimeStep(ts, dt));
  PetscCall(TSSetMaxTime(ts, .02));
  PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_INTERPOLATE));
  PetscCall(TSSetSolution(ts, x));

  if (mymonitor) {
    ctx.ports = NULL;
    PetscCall(TSMonitorSet(ts, MyMonitor, &ctx, MyDestroy));
  }

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Set runtime options
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  PetscCall(TSSetFromOptions(ts));

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Solve nonlinear system
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  PetscCall(TSSolve(ts, x));
  PetscCall(TSGetStepNumber(ts, &steps));
  PetscCall(PetscOptionsHasName(NULL, NULL, "-square_initial", &flg));
  if (flg) {
    PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "InitialSolution.heat", FILE_MODE_WRITE, &viewer));
    PetscCall(VecView(x, viewer));
    PetscCall(PetscViewerDestroy(&viewer));
  }

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Free work space.  All PETSc objects should be destroyed when they
     are no longer needed.
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  PetscCall(VecDestroy(&x));
  PetscCall(VecDestroy(&r));
  PetscCall(TSDestroy(&ts));
  PetscCall(DMDestroy(&da));

  PetscCall(PetscFinalize());
  return 0;
}
/* ------------------------------------------------------------------- */
/*
   FormFunction - Evaluates nonlinear function, F(x).

   Input Parameters:
.  ts - the TS context
.  X - input vector
.  ptr - optional user-defined context, as set by SNESSetFunction()

   Output Parameter:
.  F - function vector
 */
PetscErrorCode FormFunction(TS ts, PetscReal ftime, Vec X, Vec F, void *ptr)
{
  DM           da;
  PetscInt     i, Mx, xs, xm;
  PetscReal    hx, sx;
  PetscScalar *x, *f;
  Vec          localX;
  UserCtx     *ctx = (UserCtx *)ptr;

  PetscFunctionBegin;
  PetscCall(TSGetDM(ts, &da));
  PetscCall(DMGetLocalVector(da, &localX));
  PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));

  hx = 1.0 / (PetscReal)Mx;
  sx = 1.0 / (hx * hx);

  /*
     Scatter ghost points to local vector,using the 2-step process
        DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
     By placing code between these two statements, computations can be
     done while messages are in transition.
  */
  PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
  PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));

  /*
     Get pointers to vector data
  */
  PetscCall(DMDAVecGetArrayRead(da, localX, &x));
  PetscCall(DMDAVecGetArray(da, F, &f));

  /*
     Get local grid boundaries
  */
  PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));

  /*
     Compute function over the locally owned part of the grid
  */
  for (i = xs; i < xs + xm; i++) {
    f[i] = ctx->kappa * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx;
    if (ctx->allencahn) f[i] += (x[i] - x[i] * x[i] * x[i]);
  }

  /*
     Restore vectors
  */
  PetscCall(DMDAVecRestoreArrayRead(da, localX, &x));
  PetscCall(DMDAVecRestoreArray(da, F, &f));
  PetscCall(DMRestoreLocalVector(da, &localX));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/* ------------------------------------------------------------------- */
PetscErrorCode FormInitialSolution(DM da, Vec U)
{
  PetscInt           i, xs, xm, Mx, scale = 1, N;
  PetscScalar       *u;
  const PetscScalar *f;
  PetscReal          hx, x, r;
  Vec                finesolution;
  PetscViewer        viewer;
  PetscBool          flg;

  PetscFunctionBegin;
  PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));

  hx = 1.0 / (PetscReal)Mx;

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

  /*
     Get local grid boundaries
  */
  PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));

  /*  InitialSolution is obtained with
      ./heat -ts_monitor -snes_monitor  -pc_type lu   -snes_converged_reason    -ts_type cn  -da_refine 9 -ts_max_time 1.e-4 -ts_time_step .125e-6 -snes_atol 1.e-25 -snes_rtol 1.e-25  -ts_max_steps 15
  */
  PetscCall(PetscOptionsHasName(NULL, NULL, "-square_initial", &flg));
  if (!flg) {
    PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "InitialSolution.heat", FILE_MODE_READ, &viewer));
    PetscCall(VecCreate(PETSC_COMM_WORLD, &finesolution));
    PetscCall(VecLoad(finesolution, viewer));
    PetscCall(PetscViewerDestroy(&viewer));
    PetscCall(VecGetSize(finesolution, &N));
    scale = N / Mx;
    PetscCall(VecGetArrayRead(finesolution, &f));
  }

  /*
     Compute function over the locally owned part of the grid
  */
  for (i = xs; i < xs + xm; i++) {
    x = i * hx;
    r = PetscSqrtReal((x - .5) * (x - .5));
    if (r < .125) u[i] = 1.0;
    else u[i] = -.5;

    /* With the initial condition above the method is first order in space */
    /* this is a smooth initial condition so the method becomes second order in space */
    /*u[i] = PetscSinScalar(2*PETSC_PI*x); */
    /*  u[i] = f[scale*i];*/
    if (!flg) u[i] = f[scale * i];
  }
  if (!flg) {
    PetscCall(VecRestoreArrayRead(finesolution, &f));
    PetscCall(VecDestroy(&finesolution));
  }

  /*
     Restore vectors
  */
  PetscCall(DMDAVecRestoreArray(da, U, &u));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*
    This routine is not parallel
*/
PetscErrorCode MyMonitor(TS ts, PetscInt step, PetscReal time, Vec U, void *ptr)
{
  UserCtx            *ctx = (UserCtx *)ptr;
  PetscDrawLG         lg;
  PetscScalar        *u;
  PetscInt            Mx, i, xs, xm, cnt;
  PetscReal           x, y, hx, pause, sx, len, max, xx[2], yy[2];
  PetscDraw           draw;
  Vec                 localU;
  DM                  da;
  int                 colors[] = {PETSC_DRAW_YELLOW, PETSC_DRAW_RED, PETSC_DRAW_BLUE};
  const char *const   legend[] = {"-kappa (\\grad u,\\grad u)", "(1 - u^2)^2"};
  PetscDrawAxis       axis;
  PetscDrawViewPorts *ports;
  PetscReal           vbounds[] = {-1.1, 1.1};

  PetscFunctionBegin;
  PetscCall(PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD), 1, vbounds));
  PetscCall(PetscViewerDrawResize(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD), 1200, 800));
  PetscCall(TSGetDM(ts, &da));
  PetscCall(DMGetLocalVector(da, &localU));
  PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));
  PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
  hx = 1.0 / (PetscReal)Mx;
  sx = 1.0 / (hx * hx);
  PetscCall(DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU));
  PetscCall(DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU));
  PetscCall(DMDAVecGetArrayRead(da, localU, &u));

  PetscCall(PetscViewerDrawGetDrawLG(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD), 1, &lg));
  PetscCall(PetscDrawLGGetDraw(lg, &draw));
  PetscCall(PetscDrawCheckResizedWindow(draw));
  if (!ctx->ports) PetscCall(PetscDrawViewPortsCreateRect(draw, 1, 3, &ctx->ports));
  ports = ctx->ports;
  PetscCall(PetscDrawLGGetAxis(lg, &axis));
  PetscCall(PetscDrawLGReset(lg));

  xx[0] = 0.0;
  xx[1] = 1.0;
  cnt   = 2;
  PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-zoom", xx, &cnt, NULL));
  xs = xx[0] / hx;
  xm = (xx[1] - xx[0]) / hx;

  /*
      Plot the  energies
  */
  PetscCall(PetscDrawLGSetDimension(lg, 1 + (ctx->allencahn ? 1 : 0)));
  PetscCall(PetscDrawLGSetColors(lg, colors + 1));
  PetscCall(PetscDrawViewPortsSet(ports, 2));
  x = hx * xs;
  for (i = xs; i < xs + xm; i++) {
    xx[0] = xx[1] = x;
    yy[0]         = PetscRealPart(.25 * ctx->kappa * (u[i - 1] - u[i + 1]) * (u[i - 1] - u[i + 1]) * sx);
    if (ctx->allencahn) yy[1] = .25 * PetscRealPart((1. - u[i] * u[i]) * (1. - u[i] * u[i]));
    PetscCall(PetscDrawLGAddPoint(lg, xx, yy));
    x += hx;
  }
  PetscCall(PetscDrawGetPause(draw, &pause));
  PetscCall(PetscDrawSetPause(draw, 0.0));
  PetscCall(PetscDrawAxisSetLabels(axis, "Energy", "", ""));
  PetscCall(PetscDrawLGSetLegend(lg, legend));
  PetscCall(PetscDrawLGDraw(lg));

  /*
      Plot the  forces
  */
  PetscCall(PetscDrawViewPortsSet(ports, 1));
  PetscCall(PetscDrawLGReset(lg));
  x   = xs * hx;
  max = 0.;
  for (i = xs; i < xs + xm; i++) {
    xx[0] = xx[1] = x;
    yy[0]         = PetscRealPart(ctx->kappa * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx);
    max           = PetscMax(max, PetscAbs(yy[0]));
    if (ctx->allencahn) {
      yy[1] = PetscRealPart(u[i] - u[i] * u[i] * u[i]);
      max   = PetscMax(max, PetscAbs(yy[1]));
    }
    PetscCall(PetscDrawLGAddPoint(lg, xx, yy));
    x += hx;
  }
  PetscCall(PetscDrawAxisSetLabels(axis, "Right hand side", "", ""));
  PetscCall(PetscDrawLGSetLegend(lg, NULL));
  PetscCall(PetscDrawLGDraw(lg));

  /*
        Plot the solution
  */
  PetscCall(PetscDrawLGSetDimension(lg, 1));
  PetscCall(PetscDrawViewPortsSet(ports, 0));
  PetscCall(PetscDrawLGReset(lg));
  x = hx * xs;
  PetscCall(PetscDrawLGSetLimits(lg, x, x + (xm - 1) * hx, -1.1, 1.1));
  PetscCall(PetscDrawLGSetColors(lg, colors));
  for (i = xs; i < xs + xm; i++) {
    xx[0] = x;
    yy[0] = PetscRealPart(u[i]);
    PetscCall(PetscDrawLGAddPoint(lg, xx, yy));
    x += hx;
  }
  PetscCall(PetscDrawAxisSetLabels(axis, "Solution", "", ""));
  PetscCall(PetscDrawLGDraw(lg));

  /*
      Print the  forces as arrows on the solution
  */
  x   = hx * xs;
  cnt = xm / 60;
  cnt = (!cnt) ? 1 : cnt;

  for (i = xs; i < xs + xm; i += cnt) {
    y   = PetscRealPart(u[i]);
    len = .5 * PetscRealPart(ctx->kappa * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx) / max;
    PetscCall(PetscDrawArrow(draw, x, y, x, y + len, PETSC_DRAW_RED));
    if (ctx->allencahn) {
      len = .5 * PetscRealPart(u[i] - u[i] * u[i] * u[i]) / max;
      PetscCall(PetscDrawArrow(draw, x, y, x, y + len, PETSC_DRAW_BLUE));
    }
    x += cnt * hx;
  }
  PetscCall(DMDAVecRestoreArrayRead(da, localU, &x));
  PetscCall(DMRestoreLocalVector(da, &localU));
  PetscCall(PetscDrawStringSetSize(draw, .2, .2));
  PetscCall(PetscDrawFlush(draw));
  PetscCall(PetscDrawSetPause(draw, pause));
  PetscCall(PetscDrawPause(draw));
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscErrorCode MyDestroy(void **ptr)
{
  UserCtx *ctx = *(UserCtx **)ptr;

  PetscFunctionBegin;
  PetscCall(PetscDrawViewPortsDestroy(ctx->ports));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*TEST

   test:
     args: -ts_monitor -snes_monitor -pc_type lu -snes_converged_reason -ts_type cn -da_refine 2 -square_initial

   test:
     suffix: 2
     args: -ts_monitor -snes_monitor -pc_type lu -snes_converged_reason -ts_type cn -da_refine 5 -mymonitor -square_initial -allen-cahn -kappa .001
     requires: x

TEST*/