/***************************************
  $Header: /cvsroot/petscgraphics/cahnhill.c,v 1.8 2003/05/21 22:52:29 hazelsct Exp $

  This file contains the heart of the Cahn-Hilliard formulation, in particular
  the functions which build the finite difference residuals and Jacobian.
***************************************/


/* Including cahnhill.h includes the necessary PETSc header files */
#include "cahnhill.h"


/*++++++++++++++++++++++++++++++++++++++
  This calculates the derivative of homogeneous free energy with respect to
  concentration, which is a component of the chemical potential.  It currently
  uses
  +latex+$$\Psi' = C(1-C)\left(\frac{1}{2}+m-C\right)$$
  +html+ <center>Psi' = <i>C</i> (1-<i>C</i>) (0.5+<i>m</i>+<i>C</i>)</center>
  which gives (meta)stable equilibria at
  +latex+$C=0$ and 1 and an unstable equilibrium at $C=\frac{1}{2}+m$; if $m>0$
  +html+ <i>C</i>=0 and 1 and an unstable equilibrium at <i>C</i> = 1/2 +
  +html+ <i>m</i>; if <i>m</i>&gt;0
  then the 0 phase is stable and vice versa.

  PetscScalar ch_psiprime It returns the derivative itself.

  PetscScalar C The concentration.

  PetscScalar mparam The model parameter
  +latex+$m$.
  +html+ <i>m</i>.
  ++++++++++++++++++++++++++++++++++++++*/

static inline PetscScalar ch_psiprime (PetscScalar C, PetscScalar mparam)
{ return C*(1-C)*(0.5+mparam-C); }

#define PSIPRIME_FLOPS 5


/*++++++++++++++++++++++++++++++++++++++
  This calculates the second derivative of homogeneous free energy with respect
  to concentration, for insertion into the Jacobian matrix.  See the
  +latex+{\tt ch\_psiprime()} function for the definition of $\Psi$.
  +html+ <tt>ch_psiprime()</tt> function for the definition of Psi.

  PetscScalar ch_psidoubleprime It returns the second derivative
  +latex+$\Psi''(C)$.
  -latex-Psi''(C).

  PetscScalar C The concentration.

  PetscScalar mparam The model parameter
  +latex+$m$.
  +html+ <i>m</i>.
  ++++++++++++++++++++++++++++++++++++++*/

static inline PetscScalar ch_psidoubleprime (PetscScalar C, PetscScalar mparam)
{ return 3*C*C - (3+2*mparam)*C + 0.5+mparam; }

#define PSIDOUBLEPRIME_FLOPS 8


/*++++++++++++++++++++++++++++++++++++++
  This function computes the residual from indices to points in the
  concentration array.  ``Up'' refers to the positive
  +latex+$y$-direction, ``down'' to negative $y$, ``left'' to negative $x$ and
  +latex+``right'' to positive $x$.
  +html+ <i>y</i>-direction, ``down'' to negative <i>y</i>, ``left'' to
  +html+ negative <i>x</i> and ``right'' to positive <i>x</i>.

  PetscScalar ch_residual_2d Returns the residual itself

  PetscScalar *conc Array of concentrations

  PetscScalar alpha Model parameter
  +latex+$\kappa\alpha$
  -latex-kappa alpha

  PetscScalar beta Model parameter
  +latex+$\kappa\beta$
  -latex-kappa beta

  PetscScalar mparam Model parameter
  +latex+$m$
  -latex-m

  PetscScalar hx Inverse square
  +latex+$x$-spacing $h_x^{-2}$
  +html+ <i>x</i>-spacing <i>h<sub>x</sub></i><sup>-2</sup>

  PetscScalar hy Inverse square
  +latex+$y$-spacing $h_y^{-2}$
  +html+ <i>y</i>-spacing <i>h<sub>y</sub></i><sup>-2</sup>

  int upup Index to array position two cells up from current

  int upleft Index to array position one cell up and one left from current

  int up Index to array position one cell up from current

  int upright Index to array position one cell up and one right from current

  int leftleft Index to array position two cells left of current

  int left Index to array position one cell left of current

  int current Index to current cell array position

  int right Index to array position one cell right of current

  int rightright Index to array position two cells right of current

  int downleft Index to array position one cell down and one left from current

  int down Index to array position one cell down from current

  int downright Index to array position one cell down and one right from current

  int downdown Index to array position two cells down from current
  ++++++++++++++++++++++++++++++++++++++*/

static inline PetscScalar ch_residual_2d
(PetscScalar *conc, PetscScalar alpha, PetscScalar beta, PetscScalar mparam,
 PetscScalar hx, PetscScalar hy, int upup, int upleft, int up, int upright, int
 leftleft, int left, int current, int right, int rightright, int downleft, int
 down, int downright, int downdown)
{
  PetscScalar retval, hx2, hy2, hxhy;

  hx2 = hx*hx;  hy2 = hy*hy;  hxhy = hx*hy;

  /*+ This calculates the
    +latex+$\beta$-term, $\kappa\beta$ times the Laplacian of $\Psi'(C)$,
    -latex-beta term, kappa*beta times the laplacian of Psi prime(C),
    +*/
  retval = beta *
    (hx * (ch_psiprime(conc[left], mparam) + ch_psiprime(conc[right], mparam))
     + hy * (ch_psiprime(conc[up], mparam) + ch_psiprime(conc[down], mparam))
     - 2*(hx+hy) * ch_psiprime(conc[current], mparam));

  /*+ then subtracts the
    +latex+$\alpha$-term, $\kappa\alpha\nabla^2\nabla^2C$.
    -latex-alpha term, kappa*alpha times the Laplacian of the Laplacian of C.
    +*/
  retval -= alpha *
    (hx2 * (conc[leftleft] + conc[rightright])
     + hy2 * (conc[upup] + conc[downdown])
     + 2*hxhy * (conc[upleft]+conc[upright]+conc[downleft]+conc[downright])
     - (4*hx2 + 4*hxhy) * (conc[left] + conc[right])
     - (4*hy2 + 4*hxhy) * (conc[up] + conc[down])
     + (6*hx2 + 6*hy2 + 8*hxhy) * conc[current]);

  return retval;
}

#define RESIDUAL_FLOPS_2D (5*PSIPRIME_FLOPS + 10 + 32)


/*++++++++++++++++++++++++++++++++++++++
  This function computes the residual from indices to points in the
  concentration array.  ``Front refers to the positive
  +latex+$z$-direction, ``back'' to negative $z$, ``up'' to positive
  +latex+$y$, ``down'' to negative $y$, ``left'' to
  +latex+negative $x$ and ``right'' to positive $x$.
  +html+ <i>z</i>-direction, ``back'' to negative <i>z</i>, ``up'' to positive
  +html+ <i>y</i>, ``down'' to negative <i>y</i>, ``left'' to
  +html+ negative <i>x</i> and ``right'' to positive <i>x</i>.

  PetscScalar ch_residual_3d Returns the residual itself

  PetscScalar *conc Array of concentrations

  PetscScalar alpha Model parameter
  +latex+$\kappa\alpha$
  -latex-kappa alpha

  PetscScalar beta Model parameter
  +latex+$\kappa\beta$
  -latex-kappa beta

  PetscScalar mparam Model parameter
  +latex+$m$
  -latex-m

  PetscScalar hx Inverse square
  +latex+$x$-spacing $h_x^{-2}$
  +html+ <i>x</i>-spacing <i>h<sub>x</sub></i><sup>-2</sup>

  PetscScalar hy Inverse square
  +latex+$y$-spacing $h_y^{-2}$
  +html+ <i>y</i>-spacing <i>h<sub>y</sub></i><sup>-2</sup>

  PetscScalar hz Inverse square
  +latex+$z$-spacing $h_z^{-2}$
  +html+ <i>z</i>-spacing <i>h<sub>z</sub></i><sup>-2</sup>

  int frontfront Index to array position two cells in front of current

  int frontup Index to array position one cell front and one up from current

  int frontleft Index to array position one cell front and one left from current

  int front Index to array position one cell in front of current

  int frontright Index to array position one cell front and one right from current

  int frontdown Index to array position one cell front and one down from current

  int upup Index to array position two cells up from current

  int upleft Index to array position one cell up and one left from current

  int up Index to array position one cell up from current

  int upright Index to array position one cell up and one right from current

  int leftleft Index to array position two cells left of current

  int left Index to array position one cell left of current

  int current Index to current cell array position

  int right Index to array position one cell right of current

  int rightright Index to array position two cells right of current

  int downleft Index to array position one cell down and one left from current

  int down Index to array position one cell down from current

  int downright Index to array position one cell down and one right from current

  int downdown Index to array position two cells down from current

  int backup Index to array position one cell back and one up from current

  int backleft Index to array position one cell back and one left from current

  int back Index to array position one cell in back of current

  int backright Index to array position one cell back and one right from current

  int backdown Index to array position one cell back and one down from current

  int backback Index to array position two cells in back of current
  ++++++++++++++++++++++++++++++++++++++*/

static inline PetscScalar ch_residual_3d
(PetscScalar *conc, PetscScalar alpha, PetscScalar beta, PetscScalar mparam,
 PetscScalar hx, PetscScalar hy, PetscScalar hz, int frontfront,
 int frontup, int frontleft, int front, int frontright, int frontdown,
 int upup, int upleft, int up, int upright, int leftleft, int left,
 int current, int right, int rightright, int downleft, int down, int downright,
 int downdown,
 int backup, int backleft, int back, int backright, int backdown,
 int backback)
{
  PetscScalar retval, hx2, hy2, hz2, hxhy, hxhz, hyhz;

  hx2 = hx*hx;  hy2 = hy*hy;  hz2 = hz*hz;
  hxhy = hx*hy; hxhz = hx*hz; hyhz = hy*hz;

  /*+ This calculates the
    +latex+$\beta$-term, $\kappa\beta$ times the Laplacian of $\Psi'(C)$,
    -latex-beta term, kappa*beta times the laplacian of Psi prime(C),
    +*/
  retval = beta *
    (hx * (ch_psiprime(conc[left], mparam) + ch_psiprime(conc[right], mparam))
     + hy * (ch_psiprime(conc[up], mparam) + ch_psiprime(conc[down], mparam))
     + hz * (ch_psiprime(conc[front], mparam)+ ch_psiprime(conc[back], mparam))
     - 2*(hx+hy+hz) * ch_psiprime(conc[current], mparam));

  /*+ then subtracts the
    +latex+$\alpha$-term, $\kappa\alpha\nabla^2\nabla^2C$.
    -latex-alpha term, kappa*alpha times the Laplacian of the Laplacian of C.
    +*/
  retval -= alpha *
    (hx2 * (conc[leftleft] + conc[rightright])
     + hy2 * (conc[upup] + conc[downdown])
     + hz2 * (conc[frontfront] + conc[backback])
     + 2*hxhy * (conc[upleft]+conc[upright]+conc[downleft]+conc[downright])
     + 2*hxhz*(conc[frontleft]+conc[frontright]+conc[backleft]+conc[backright])
     + 2*hyhz * (conc[frontup]+conc[frontdown]+conc[backup]+conc[backdown])
     - (4*hx2 + 4*hxhy + 4*hxhz) * (conc[left] + conc[right])
     - (4*hy2 + 4*hxhy + 4*hyhz) * (conc[up] + conc[down])
     - (4*hz2 + 4*hxhz + 4*hyhz) * (conc[front] + conc[back])
     + (6*hx2 + 6*hy2 + 6*hz2 + 8*hxhy + 8*hxhz + 8*hyhz) * conc[current]);

  return retval;
}

#define RESIDUAL_FLOPS_3D (7*PSIPRIME_FLOPS + 14 + 65)


/*++++++++++++++++++++++++++++++++++++++
  This evaluates the nonlinear finite difference approximation to the residuals
  +latex+$R_i$.
  +html+ <i>R<sub>i</sub></i>.
  Note that it loops on the interior points and the boundary separately, to
  avoid conditional statements within the double loop over the local grid
  indices.

  int ch_residual_vector_2d It returns zero or an error value.

  Vec X The pre-allocated local vector of unknowns.

  Vec F The pre-allocated local vector of residuals, filled by this function.

  void *ptr Data passed in the application context.
  ++++++++++++++++++++++++++++++++++++++*/

int ch_residual_vector_2d (Vec X, Vec F, void *ptr)
{
  AppCtx  *user = (AppCtx*)ptr;
  int     ierr,i,j, mc,chvar, mx,my,xs,ys,xm,ym,gxs,gys,gxm,gym;
  int     xints,xinte,yints,yinte;
  PetscScalar  hx,hy,dhx,dhy,hxm2,hym2;
  PetscScalar  alpha, beta, mparam;
  PetscScalar  *x,*f;
  Vec     localX = user->localX,localF = user->localF; 

  /* Scatter ghost points to local vector, using the 2-step process
        DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
     By placing code between these two statements, computations can be
     done while messages are in transition. */
  ierr = DAGlobalToLocalBegin (user->da, X, INSERT_VALUES, localX);
  CHKERRQ (ierr);

  /* Define mesh intervals ratios for uniform grid.
     [Note: FD formulae below may someday be normalized by multiplying through
     by local volume element to obtain coefficients O(1) in two dimensions.] */

  mc = user->mc;
  chvar = user->chvar;
  mx = user->mx;
  my = user->my;
  mparam = user->mparam;
  alpha = user->kappa * user->gamma * user->epsilon;
  beta = user->kappa * user->gamma / user->epsilon;
  dhx = (PetscScalar)(mx-1);
  dhy = (PetscScalar)(my-1);
  /* These next two lines hard-code the 1x1 square */
  hx = 1./dhx;
  hy = 1./dhy;
  hxm2 = 1./hx/hx;
  hym2 = 1./hy/hy;

  ierr = DAGlobalToLocalEnd (user->da, X, INSERT_VALUES, localX);
  CHKERRQ (ierr);

  /* Get pointers to vector data */
  ierr = VecGetArray (localX, &x); CHKERRQ (ierr);
  ierr = VecGetArray (localF, &f); CHKERRQ (ierr);

  /* Get local grid boundaries */
  ierr = DAGetCorners (user->da, &xs, &ys, PETSC_NULL, &xm, &ym, PETSC_NULL);
  CHKERRQ (ierr);
  ierr = DAGetGhostCorners (user->da, &gxs, &gys, PETSC_NULL, &gxm, &gym,
			    PETSC_NULL); CHKERRQ (ierr);

  /* Determine the interior of the locally owned part of the grid. */
  if (user->period == DA_XPERIODIC || user->period == DA_XYPERIODIC) {
    xints = xs; xinte = xs+xm; }
  else {
    xints = (xs==0) ? 2:xs; xinte = (xs+xm==mx) ? xs+xm-2 : xs+xm; }
  if (user->period == DA_YPERIODIC || user->period == DA_XYPERIODIC) {
    yints = ys; yinte = ys+ym; }
  else {
    yints = (ys==0) ? 2:ys; yinte = (ys+ym==my) ? ys+ym-2 : ys+ym; }

  /* Loop over the interior points */
  for (j=yints-gys; j<yinte-gys; j++)
    for (i=j*gxm + xints-gxs; i<j*gxm + xinte-gxs; i++)
      f[C(i)] = ch_residual_2d
	(x, alpha, beta, mparam, hxm2, hym2,
	 C(i+2*gxm), C(i+gxm-1), C(i+gxm), C(i+gxm+1),
	 C(i-2), C(i-1), C(i), C(i+1), C(i+2),
	 C(i-gxm-1), C(i-gxm), C(i-gxm+1), C(i-2*gxm));

  /* Okay, that was the easy part.  Now for the fun part.
     The following conditionals implement symmetry boundary conditions. */

  /* Test whether we are on the bottom edge of the global array */
  if (yints == 2) /* Not ys==0 because that would be true for periodic too */
    {
      printf ("This must only be reached if we are NOT y-periodic\n");
      /* Bottom two edge lines */
      for (i=xints-gxs; i<xinte-gxs; i++)
	{
	  f[C(i)] = ch_residual_2d
	    (x, alpha, beta, mparam, hxm2, hym2,
	     C(i+2*gxm), C(i+gxm-1), C(i+gxm), C(i+gxm+1),
	     C(i-2), C(i-1), C(i), C(i+1), C(i+2),
	     C(i+gxm-1), C(i+gxm), C(i+gxm+1), C(i+2*gxm));
	  f[C(i+gxm)] = ch_residual_2d
	    (x, alpha, beta, mparam, hxm2, hym2,
	     C(i+3*gxm), C(i+2*gxm-1), C(i+2*gxm), C(i+2*gxm+1),
	     C(i+gxm-2), C(i+gxm-1), C(i+gxm), C(i+gxm+1), C(i+gxm+2),
	     C(i-1), C(i), C(i+1), C(i+gxm));
	}

      /* Bottom left corner */
      if (xints == 2)
	{
	  f[C(0)] = ch_residual_2d
	    (x, alpha, beta, mparam, hxm2, hym2,
	     C(2*gxm), C(gxm+1), C(gxm), C(gxm+1),
	     C(2), C(1), C(0), C(1), C(2),
	     C(gxm+1), C(gxm), C(gxm+1), C(2*gxm));
	  f[C(1)] = ch_residual_2d
	    (x, alpha, beta, mparam, hxm2, hym2,
	     C(2*gxm+1), C(gxm), C(gxm+1), C(gxm+2),
	     C(1), C(0), C(1), C(2), C(3),
	     C(gxm), C(gxm+1), C(gxm+2), C(2*gxm+1));
	  f[C(gxm)] = ch_residual_2d
	    (x, alpha, beta, mparam, hxm2, hym2,
	     C(3*gxm), C(2*gxm+1), C(2*gxm), C(2*gxm+1),
	     C(gxm+2), C(gxm+1), C(gxm+0), C(gxm+1), C(gxm+2),
	     C(1), C(0), C(1), C(gxm));
	  f[C(gxm+1)] = ch_residual_2d
	    (x, alpha, beta, mparam, hxm2, hym2,
	     C(3*gxm+1), C(2*gxm), C(2*gxm+1), C(2*gxm+2),
	     C(gxm+1), C(gxm), C(gxm+1), C(gxm+2), C(gxm+3),
	     C(0), C(1), C(2), C(gxm+1));
	}

      /* Bottom right corner */
      if (xinte == mx-2)
	{
	  i = mx-gxs-1; /* Array index of the bottom right corner point */
	  f[C(i)] = ch_residual_2d
	    (x, alpha, beta, mparam, hxm2, hym2,
	     C(i+2*gxm), C(i+gxm-1), C(i+gxm), C(i+gxm-1),
	     C(i-2), C(i-1), C(i), C(i-1), C(i-2),
	     C(i+gxm-1), C(i+gxm), C(i+gxm-1), C(i+2*gxm));
	  f[C(i-1)] = ch_residual_2d
	    (x, alpha, beta, mparam, hxm2, hym2,
	     C(i+2*gxm-1), C(i+gxm-2), C(i+gxm-1), C(i+gxm),
	     C(i-3), C(i-2), C(i-1), C(i), C(i-1),
	     C(i+gxm-2), C(i+gxm-1), C(i+gxm), C(i+2*gxm-1));
	  f[C(i+gxm)] = ch_residual_2d
	    (x, alpha, beta, mparam, hxm2, hym2,
	     C(i+3*gxm), C(i+2*gxm-1), C(i+2*gxm), C(i+2*gxm-1),
	     C(i+gxm-2), C(i+gxm-1), C(i+gxm), C(i+gxm-1), C(i+gxm-2),
	     C(i-1), C(i), C(i-1), C(i+gxm));
	  f[C(i+gxm-1)] = ch_residual_2d
	    (x, alpha, beta, mparam, hxm2, hym2,
	     C(i+3*gxm-1), C(i+2*gxm-2), C(i+2*gxm-1), C(i+2*gxm),
	     C(i+gxm-3), C(i+gxm-2), C(i+gxm-1), C(i+gxm), C(i+gxm-1),
	     C(i-2), C(i-1), C(i), C(i+gxm-1));
	}
    }

  /* Test whether we are on the top edge of the global array */
  if (yinte == my-2)
    {
      printf ("This must only be reached if we are NOT y-periodic\n");
      /* Top two edge lines */
      for (i=(my-gys-1)*gxm + xints-gxs; i<(my-gys-1)*gxm + xinte-gxs; i++)
	{
	  f[C(i)] = ch_residual_2d
	    (x, alpha, beta, mparam, hxm2, hym2,
	     C(i-2*gxm), C(i-gxm-1), C(i-gxm), C(i-gxm+1),
	     C(i-2), C(i-1), C(i), C(i+1), C(i+2),
	     C(i-gxm-1), C(i-gxm), C(i-gxm+1), C(i-2*gxm));
	  f[C(i-gxm)] = ch_residual_2d
	    (x, alpha, beta, mparam, hxm2, hym2,
	     C(i-gxm), C(i-1), C(i), C(i+1),
	     C(i-gxm-2), C(i-gxm-1), C(i-gxm), C(i-gxm+1), C(i-gxm+2),
	     C(i-2*gxm-1), C(i-2*gxm), C(i-2*gxm+1), C(i-3*gxm));
	}

      /* Top left corner */
      if (xints == 2)
	{
	  i = (my-gys-1)*gxm; /* Array index of the top left corner point */
	  f[C(i)] = ch_residual_2d
	    (x, alpha, beta, mparam, hxm2, hym2,
	     C(i-2*gxm), C(i-gxm+1), C(i-gxm), C(i-gxm+1),
	     C(i+2), C(i+1), C(i), C(i+1), C(i+2),
	     C(i-gxm+1), C(i-gxm), C(i-gxm+1), C(i-2*gxm));
	  f[C(i+1)] = ch_residual_2d
	    (x, alpha, beta, mparam, hxm2, hym2,
	     C(i-2*gxm+1), C(i-gxm), C(i-gxm+1), C(i-gxm+2),
	     C(i+1), C(i), C(i+1), C(i+2), C(i+3),
	     C(i-gxm), C(i-gxm+1), C(i-gxm+2), C(i-2*gxm+1));
	  f[C(i-gxm)] = ch_residual_2d
	    (x, alpha, beta, mparam, hxm2, hym2,
	     C(i-gxm), C(i+1), C(i), C(i+1),
	     C(i-gxm+2), C(i-gxm+1), C(i-gxm), C(i-gxm+1), C(i-gxm+2),
	     C(i-2*gxm+1), C(i-2*gxm), C(i-2*gxm+1), C(i-3*gxm));
	  f[C(i-gxm+1)] = ch_residual_2d
	    (x, alpha, beta, mparam, hxm2, hym2,
	     C(i-gxm+1), C(i), C(i+1), C(i+2),
	     C(i-gxm+1), C(i-gxm), C(i-gxm+1), C(i-gxm+2), C(i-gxm+3),
	     C(i-2*gxm), C(i-2*gxm+1), C(i-2*gxm+2), C(i-3*gxm+1));
	}

      /* Top right corner */
      if (xinte == mx-2)
	{ /* Array index of top right corner point */
	  i = (my-gys-1)*gxm + mx-gxs-1;
	  f[C(i)] = ch_residual_2d
	    (x, alpha, beta, mparam, hxm2, hym2,
	     C(i-2*gxm), C(i-gxm-1), C(i-gxm), C(i-gxm-1),
	     C(i-2), C(i-1), C(i), C(i-1), C(i-2),
	     C(i-gxm-1), C(i-gxm), C(i-gxm-1), C(i-2*gxm));
	  f[C(i-1)] = ch_residual_2d
	    (x, alpha, beta, mparam, hxm2, hym2,
	     C(i-2*gxm-1), C(i-gxm-2), C(i-gxm-1), C(i-gxm),
	     C(i-3), C(i-2), C(i-1), C(i), C(i-1),
	     C(i-gxm-2), C(i-gxm-1), C(i-gxm), C(i-2*gxm-1));
	  f[C(i-gxm)] = ch_residual_2d
	    (x, alpha, beta, mparam, hxm2, hym2,
	     C(i-gxm), C(i-1), C(i), C(i-1),
	     C(i-gxm-2), C(i-gxm-1), C(i-gxm), C(i-gxm-1), C(i-gxm-2),
	     C(i-2*gxm-1), C(i-2*gxm), C(i-2*gxm-1), C(i-3*gxm));
	  f[C(i-gxm-1)] = ch_residual_2d
	    (x, alpha, beta, mparam, hxm2, hym2,
	     C(i-gxm-1), C(i-2), C(i-1), C(i),
	     C(i-gxm-3), C(i-gxm-2), C(i-gxm-1), C(i-gxm), C(i-gxm-1),
	     C(i-2*gxm-2), C(i-2*gxm-1), C(i-2*gxm), C(i-3*gxm-1));
	}
    }

  /* Test whether we are on the left edge of the global array */
  if (xints == 2)
    {
      printf ("This must only be reached if we are NOT x-periodic\n");
      /* Left 2 edge lines */
      for (j=yints-gys; j<yinte-gys; j++)
	{
	  i = j*gxm;
	  f[C(i)] = ch_residual_2d
	    (x, alpha, beta, mparam, hxm2, hym2,
	     C(i+2*gxm), C(i+gxm+1), C(i+gxm), C(i+gxm+1),
	     C(i+2), C(i+1), C(i), C(i+1), C(i+2),
	     C(i-gxm+1), C(i-gxm), C(i-gxm+1), C(i-2*gxm));
	  f[C(i+1)] = ch_residual_2d
	    (x, alpha, beta, mparam, hxm2, hym2,
	     C(i+2*gxm+1), C(i+gxm), C(i+gxm+1), C(i+gxm+2),
	     C(i+1), C(i), C(i+1), C(i+2), C(i+3),
	     C(i-gxm), C(i-gxm+1), C(i-gxm+2), C(i-2*gxm+1));
	}
    }

  /* Test whether we are on the right edge of the global array */
  if (xinte == mx-2)
    {
      printf ("This must only be reached if we are NOT x-periodic\n");
      /* Right 2 edge lines */ 
      for (j=yints-gys; j<yinte-gys; j++)
	{
	  i = j*gxm + mx-gxs-1;
	  f[C(i)] = ch_residual_2d
	    (x, alpha, beta, mparam, hxm2, hym2,
	     C(i+2*gxm), C(i+gxm-1), C(i+gxm), C(i+gxm-1),
	     C(i-2), C(i-1), C(i), C(i-1), C(i-2),
	     C(i-gxm-1), C(i-gxm), C(i-gxm-1), C(i-2*gxm));
	  f[C(i-1)] = ch_residual_2d
	    (x, alpha, beta, mparam, hxm2, hym2,
	     C(i+2*gxm-1), C(i+gxm-2), C(i+gxm-1), C(i+gxm),
	     C(i-3), C(i-2), C(i-1), C(i), C(i-1),
	     C(i-gxm-2), C(i-gxm-1), C(i-gxm), C(i-2*gxm-1));
	}
    }

  /* Restore vectors */
  ierr = VecRestoreArray (localX, &x); CHKERRQ (ierr);
  ierr = VecRestoreArray (localF, &f); CHKERRQ (ierr);

  /* Insert values into global vector */
  ierr = DALocalToGlobal(user->da,localF,INSERT_VALUES,F); CHKERRQ (ierr);

  /* Flop count (multiply-adds are counted as 2 operations) */
  ierr = PetscLogFlops(RESIDUAL_FLOPS_2D*ym*xm); CHKERRQ (ierr);

  return 0; 
}


/*++++++++++++++++++++++++++++++++++++++
  This evaluates the nonlinear finite difference approximation to the residuals
  +latex+$R_i$.
  +html+ <i>R<sub>i</sub></i>.
  Note that it loops on the interior points and the boundary separately, to
  avoid conditional statements within the double loop over the local grid
  indices.
  +latex+\par
  +html+ <p>
  Thus far, only periodic boundary conditions work.

  int ch_residual_vector_3d It returns zero or an error value.

  Vec X The pre-allocated local vector of unknowns.

  Vec F The pre-allocated local vector of residuals, filled by this function.

  void *ptr Data passed in the application context.
  ++++++++++++++++++++++++++++++++++++++*/

int ch_residual_vector_3d (Vec X, Vec F, void *ptr)
{
  AppCtx  *user = (AppCtx*)ptr;
  int     ierr,i,j,k, mc,chvar, mx,my,mz,xs,ys,zs,xm,ym,zm;
  int     gxs,gys,gzs,gxm,gym,gzm, xints,xinte,yints,yinte,zints,zinte;
  PetscScalar  hx,hy,hz,dhx,dhy,dhz,hxm2,hym2,hzm2;
  PetscScalar  alpha, beta, mparam;
  PetscScalar  *x,*f;
  Vec     localX = user->localX,localF = user->localF; 

  /* Scatter ghost points to local vector, using the 2-step process
        DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
     By placing code between these two statements, computations can be
     done while messages are in transition. */
  ierr = DAGlobalToLocalBegin (user->da, X, INSERT_VALUES, localX);
  CHKERRQ (ierr);

  /* Define mesh intervals ratios for uniform grid.
     [Note: FD formulae below may someday be normalized by multiplying through
     by local volume element to obtain coefficients O(1) in two dimensions.] */

  mc = user->mc;
  chvar = user->chvar;
  mx = user->mx; my = user->my; mz = user->mz;
  mparam = user->mparam;
  alpha = user->kappa * user->gamma * user->epsilon;
  beta = user->kappa * user->gamma / user->epsilon;
  dhx = (PetscScalar)(mx-1); dhy = (PetscScalar)(my-1);
  dhz = (PetscScalar)(mz-1);
  /* This next line hard-codes the 1x1x1 cube */
  hx = 1./dhx;     hy = 1./dhy;     hz = 1./dhz;
  hxm2 = 1./hx/hx; hym2 = 1./hy/hy; hzm2 = 1./hz/hz;

  ierr = DAGlobalToLocalEnd (user->da, X, INSERT_VALUES, localX);
  CHKERRQ (ierr);

  /* Get pointers to vector data */
  ierr = VecGetArray (localX, &x); CHKERRQ (ierr);
  ierr = VecGetArray (localF, &f); CHKERRQ (ierr);

  /* Get local grid boundaries */
  ierr = DAGetCorners (user->da, &xs, &ys, &zs, &xm, &ym, &zm); CHKERRQ (ierr);
  ierr = DAGetGhostCorners (user->da, &gxs, &gys, &gzs, &gxm, &gym, &gzm);
  CHKERRQ (ierr);

  /* Determine the interior of the locally owned part of the grid. */
  if (user->period == DA_XYZPERIODIC) {
    xints = xs; xinte = xs+xm;
    yints = ys; yinte = ys+ym;
    zints = zs; zinte = zs+zm; }
  else {
    SETERRQ(PETSC_ERR_ARG_WRONGSTATE,
	    "Only periodic boundary conditions in 3-D Cahn-Hilliard"); }

  /* Loop over the interior points (no boundaries yet) */
  for (k=zints-gzs; k<zinte-gzs; k++)
    for (j=k*gym + yints-gys; j<k*gym + yinte-gys; j++)
      for (i=j*gxm + xints-gxs; i<j*gxm + xinte-gxs; i++)
	f[C(i)] = ch_residual_3d
	  (x, alpha, beta, mparam, hxm2, hym2, hzm2,
	   C(i+2*gxm*gym), C(i+gxm*gym+gxm), C(i+gxm*gym-1), C(i+gxm*gym),
	   C(i+gxm*gym+1), C(i+gxm*gym-gxm),
	   C(i+2*gxm), C(i+gxm-1), C(i+gxm), C(i+gxm+1),
	   C(i-2), C(i-1), C(i), C(i+1), C(i+2),
	   C(i-gxm-1), C(i-gxm), C(i-gxm+1), C(i-2*gxm),
	   C(i-gxm*gym+gxm), C(i-gxm*gym-1), C(i-gxm*gym), C(i-gxm*gym+1),
	   C(i-gxm*gym-gxm), C(i-2*gxm*gym));

  /* Restore vectors */
  ierr = VecRestoreArray (localX, &x); CHKERRQ (ierr);
  ierr = VecRestoreArray (localF, &f); CHKERRQ (ierr);

  /* Insert values into global vector */
  ierr = DALocalToGlobal(user->da,localF,INSERT_VALUES,F); CHKERRQ (ierr);

  /* Flop count (multiply-adds are counted as 2 operations) */
  ierr = PetscLogFlops(RESIDUAL_FLOPS_3D*ym*xm*zm); CHKERRQ (ierr);

  /* ierr = VecView (F, VIEWER_STDOUT_SELF); CHKERRQ (ierr); */

  return 0;
}


/*++++++++++++++++++++++++++++++++++++++
  This computes the Jacobian matrix at each iteration, starting with the alpha
  term which is pre-computed at the beginning by
  +latex+{\tt ch\_jacobian\_alpha\_2d()}.
  +html+ <tt>ch_jacobian_alpha_2d()</tt>.

  int ch_jacobian_2d It returns 0 or an error code.

  Vec X The vector of unknowns.

  Mat *A The Jacobian matrix returned to PETSc.

  Mat *B The matrix preconditioner, in this case the same matrix.

  MatStructure *flag Flag saying the nonzeroes are the same each time.

  void *ptr Application context structure.
  ++++++++++++++++++++++++++++++++++++++*/

int ch_jacobian_2d (Vec X, Mat *A, Mat *B, MatStructure *flag, void *ptr)
{
  AppCtx  *user = (AppCtx*)ptr;
  int     ierr,node,i,j, mc,chvar, mx,my,xs,ys,xm,ym,gxs,gys,gxm,gym;
  int     xints,xinte,yints,yinte;
  int     columns [5];
  PetscScalar  hx,hy,dhx,dhy,hxm2,hym2;
  PetscScalar  alpha,beta,mparam;
  PetscScalar  *x, bvalue[5];
  Vec     localX = user->localX;

  /* Set the MatStructure flag right, Mats to return */
  *flag = SAME_NONZERO_PATTERN;
  A = &(user->J);
  B = &(user->J);

  /* Copy the alpha term Jacobian into the main Jacobian space */
  ierr = MatCopy (user->alpha, user->J, SAME_NONZERO_PATTERN);

  /* Scatter ghost points to local vector, using 2-step async I/O with (a
     couple of trivial) computations in between. */

  ierr = DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX); CHKERRQ (ierr);
  mc = user->mc;  chvar = user->chvar;
  mx = user->mx;  my = user->my;  mparam = user->mparam;
  alpha = user->kappa * user->gamma * user->epsilon;
  beta = user->kappa * user->gamma / user->epsilon;

  /* Define mesh intervals ratios for uniform grid.
     [Note: FD formulae below may someday be normalized by multiplying through
     by local volume element to obtain coefficients O(1) in two dimensions.] */
  dhx = (PetscScalar)(mx-1);  dhy = (PetscScalar)(my-1);
  hx = 1./dhx;           hy = 1./dhy; /* This line hard-codes the 1x1 square */
  hxm2 = 1./hx/hx;       hym2 = 1./hy/hy;

  ierr = DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX); CHKERRQ (ierr);

  /* Get pointer to vector data */
  ierr = VecGetArray(localX,&x); CHKERRQ (ierr);

  /*
     Get local grid boundaries
  */
  ierr = DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ (ierr);
  ierr = DAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); CHKERRQ (ierr);

  /* Determine the interior of the locally owned part of the grid. */
  if (user->period == DA_XPERIODIC || user->period == DA_XYPERIODIC) {
    xints = xs; xinte = xs+xm; }
  else {
    xints = (xs==0) ? 1:xs; xinte = (xs+xm==mx) ? xs+xm-1 : xs+xm; }
  if (user->period == DA_YPERIODIC || user->period == DA_XYPERIODIC) {
    yints = ys; yinte = ys+ym; }
  else {
    yints = (ys==0) ? 1:ys; yinte = (ys+ym==my) ? ys+ym-1 : ys+ym; }

  /* Loop over the interior points... */
  for (j=yints-gys; j<yinte-gys; j++)
    for (i=j*gxm + xints-gxs; i<j*gxm + xinte-gxs; i++) {

      /* Form the column indices */
      columns[0]=C(i+gxm);
      columns[1]=C(i-1); columns[2]=C(i); columns[3]=C(i+1);
      columns[4]=C(i-gxm);

      bvalue[0] = beta * hym2 * ch_psidoubleprime (x[columns[0]], mparam);
      bvalue[1] = beta * hxm2 * ch_psidoubleprime (x[columns[1]], mparam);
      bvalue[2] = -2.*beta * (hxm2+hym2) *
	ch_psidoubleprime (x[columns[2]], mparam);
      bvalue[3] = beta * hxm2 * ch_psidoubleprime (x[columns[3]], mparam);
      bvalue[4] = beta * hym2 * ch_psidoubleprime (x[columns[4]], mparam);

      node = C(i);
      MatSetValuesLocal (user->J, 1,&node, 5,columns, bvalue, ADD_VALUES);
    }

  /* Now the boundary conditions... */

  /* Test whether we're on the bottom row and non-y-periodic */
  if (yints == 1) {
    printf ("This must only be reached if we are NOT y-periodic\n");
    for (i=xints-gxs; i<xinte-gxs; i++) {

      /* Form the column indices */
      columns[0]=C(i+gxm);
      columns[1]=C(i-1); columns[2]=C(i); columns[3]=C(i+1);

      bvalue[0] = 2.*beta * hym2 * ch_psidoubleprime (x[columns[0]], mparam);
      bvalue[1] = beta * hxm2 * ch_psidoubleprime (x[columns[1]], mparam);
      bvalue[2] = -2.*beta * (hxm2+hym2) *
	ch_psidoubleprime (x[columns[2]], mparam);
      bvalue[3] = beta * hxm2 * ch_psidoubleprime (x[columns[3]], mparam);

      node = C(i);
      MatSetValuesLocal (user->J, 1,&node, 4,columns, bvalue, ADD_VALUES);
    }

    /* Bottom left corner */
    if (xints == 1) {
      columns[0]=C(i=0); columns[1]=C(i+1);
      columns[2]=C(i+gxm);

      bvalue[0] = -2.*beta * (hxm2+hym2) *
	ch_psidoubleprime (x[columns[0]], mparam);
      bvalue[1] = 2.*beta * hxm2 * ch_psidoubleprime (x[columns[1]], mparam);
      bvalue[2] = 2.*beta * hym2 * ch_psidoubleprime (x[columns[2]], mparam);

      node = C(i);
      MatSetValuesLocal (user->J, 1,&node, 3,columns, bvalue, ADD_VALUES);
    }

    /* Bottom right corner */
    if (xinte == mx-1) {
      columns[0]=C(i=mx-gxs-1); columns[1]=C(i-1);
      columns[2]=C(i+gxm);

      bvalue[0] = -2.*beta * (hxm2+hym2) *
	ch_psidoubleprime (x[columns[0]], mparam);
      bvalue[1] = 2.*beta * hxm2 * ch_psidoubleprime (x[columns[1]], mparam);
      bvalue[2] = 2.*beta * hym2 * ch_psidoubleprime (x[columns[2]], mparam);

      node = C(i);
      MatSetValuesLocal (user->J, 1,&node, 3,columns, bvalue, ADD_VALUES);
    }
  }

  /* Test whether we're on the top row and non-y-periodic */
  if (yinte == my-1) {
    printf ("This must only be reached if we are NOT y-periodic\n");
    for (i=(my-gys-1)*gxm + xints-gxs; i<(my-gys-1)*gxm + xinte-gxs; i++) {

      /* Form the column indices */
      columns[0]=C(i-1); columns[1]=C(i); columns[2]=C(i+1);
      columns[3]=C(i-gxm);

      bvalue[0] = beta * hxm2 * ch_psidoubleprime (x[columns[0]], mparam);
      bvalue[1] = -2.*beta * (hxm2+hym2) *
	ch_psidoubleprime (x[columns[1]], mparam);
      bvalue[2] = beta * hxm2 * ch_psidoubleprime (x[columns[2]], mparam);
      bvalue[3] = 2.*beta * hym2 * ch_psidoubleprime (x[columns[3]], mparam);

      node = C(i);
      MatSetValuesLocal (user->J, 1,&node, 4,columns, bvalue, ADD_VALUES);
    }

    /* Top left corner */
    if (xints == 1) {
      columns[0]=C(i=(my-gys-1)*gxm); columns[1] = C(i+1);
      columns[2]=C(i-gxm);

      bvalue[0] = -2.*beta * (hxm2+hym2) *
	ch_psidoubleprime (x[columns[0]], mparam);
      bvalue[1] = 2.*beta * hxm2 * ch_psidoubleprime (x[columns[1]], mparam);
      bvalue[2] = 2.*beta * hym2 * ch_psidoubleprime (x[columns[2]], mparam);

      node = C(i);
      MatSetValuesLocal (user->J, 1,&node, 3,columns, bvalue, ADD_VALUES);
    }

    /* Top right corner */
    if (xinte == mx-1) {
      columns[0]=C(i=(my-gys-1)*gxm + mx-gxs-1); columns[1]=C(i-1);
      columns[2]=C(i-gxm);

      bvalue[0] = -2.*beta * (hxm2+hym2) *
	ch_psidoubleprime (x[columns[0]], mparam);
      bvalue[1] = 2.*beta * hxm2 * ch_psidoubleprime (x[columns[1]], mparam);
      bvalue[2] = 2.*beta * hym2 * ch_psidoubleprime (x[columns[2]], mparam);

      node = C(i);
      MatSetValuesLocal (user->J, 1,&node, 3,columns, bvalue, ADD_VALUES);
    }
  }

  /* Left column */
  if (xints == 1) {
    printf ("This must only be reached if we are NOT x-periodic\n");
    for (i=(yints-gys)*gxm; i<(yinte-gys)*gxm; i+=gxm) {
      columns[0]=C(i+gxm);
      columns[1]=C(i); columns[2] = C(i+1);
      columns[3]=C(i-gxm);

      bvalue[0] = beta * hym2 * ch_psidoubleprime (x[columns[0]], mparam);
      bvalue[1] = -2.*beta * (hxm2+hym2) *
	ch_psidoubleprime (x[columns[1]], mparam);
      bvalue[2] = 2.*beta * hxm2 * ch_psidoubleprime (x[columns[2]], mparam);
      bvalue[3] = beta * hym2 * ch_psidoubleprime (x[columns[3]], mparam);

      node = C(i);
      MatSetValuesLocal (user->J, 1,&node, 4,columns, bvalue, ADD_VALUES);
    }
  }

  /* Right column */
  if (xinte == mx-1) {
    printf ("This must only be reached if we are NOT x-periodic\n");
    for (i=(yints-gys)*gxm + mx-gxs-1; i<(yinte-gys)*gxm + mx-gxs-1; i+=gxm) {
      columns[0]=C(i+gxm);
      columns[1]=C(i-1); columns[2] = C(i);
      columns[3]=C(i-gxm);

      bvalue[0] = beta * hym2 * ch_psidoubleprime (x[columns[0]], mparam);
      bvalue[1] = 2.*beta * hxm2 * ch_psidoubleprime (x[columns[1]], mparam);
      bvalue[2] = -2.*beta * (hxm2+hym2) *
	ch_psidoubleprime (x[columns[2]], mparam);
      bvalue[3] = beta * hym2 * ch_psidoubleprime (x[columns[3]], mparam);

      node = C(i);
      MatSetValuesLocal (user->J, 1,&node, 4,columns, bvalue, ADD_VALUES);
    }
  }

  /* Assemble matrix, using the 2-step process:
     MatAssemblyBegin(), MatAssemblyEnd(). */
  ierr = MatAssemblyBegin(user->J,MAT_FINAL_ASSEMBLY); CHKERRQ (ierr);
  ierr = MatAssemblyEnd(user->J,MAT_FINAL_ASSEMBLY); CHKERRQ (ierr);

  /* Restore unknown vector */
  ierr = VecRestoreArray(localX,&x); CHKERRQ (ierr);

  return ierr;
}


/*++++++++++++++++++++++++++++++++++++++
  This computes the Jacobian matrix at each iteration, starting with the alpha
  term which is pre-computed at the beginning by
  +latex+{\tt ch\_jacobian\_alpha\_3d()}.
  +html+ <tt>ch_jacobian_alpha_3d()</tt>.

  int ch_jacobian_3d It returns 0 or an error code.

  Vec X The vector of unknowns.

  Mat *A The Jacobian matrix returned to PETSc.

  Mat *B The matrix preconditioner, in this case the same matrix.

  MatStructure *flag Flag saying the nonzeroes are the same each time.

  void *ptr Application context structure.
  ++++++++++++++++++++++++++++++++++++++*/

int ch_jacobian_3d (Vec X, Mat *A, Mat *B, MatStructure *flag, void *ptr)
{
  AppCtx  *user = (AppCtx*)ptr;
  int     ierr,node,i,j,k, mc,chvar, mx,my,mz, xs,ys,zs, xm,ym,zm;
  int     gxs,gys,gzs, gxm,gym,gzm, xints,xinte, yints,yinte, zints,zinte;
  int     columns [7];
  PetscScalar  hx,hy,hz, dhx,dhy,dhz, hxm2,hym2,hzm2;
  PetscScalar  alpha,beta,mparam;
  PetscScalar  *x, bvalue[7];
  Vec     localX = user->localX;

  /* Set the MatStructure flag right, Mats to return */
  *flag = SAME_NONZERO_PATTERN;
  A = &(user->J);
  B = &(user->J);

  /* Copy the alpha term Jacobian into the main Jacobian space */
  ierr = MatCopy (user->alpha, user->J, SAME_NONZERO_PATTERN);

  /* Scatter ghost points to local vector, using 2-step async I/O with (a
     couple of trivial) computations in between. */

  ierr = DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX); CHKERRQ (ierr);
  mc = user->mc;  chvar = user->chvar;
  mx = user->mx;  my = user->my;  mz = user->mz;  mparam = user->mparam;
  alpha = user->kappa * user->gamma * user->epsilon;
  beta = user->kappa * user->gamma / user->epsilon;

  /* Define mesh intervals ratios for uniform grid.
     [Note: FD formulae below may someday be normalized by multiplying through
     by local volume element to obtain coefficients O(1) in two dimensions.] */
  dhx = (PetscScalar)(mx-1);  dhy = (PetscScalar)(my-1);
  dhz = (PetscScalar)(mz-1);
  /* This line hard-codes the 1x1x1 cube */
  hx = 1./dhx;           hy = 1./dhy;          hz = 1./dhz;
  hxm2 = 1./hx/hx;       hym2 = 1./hy/hy;      hzm2 = 1./hz/hz;

  ierr = DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX); CHKERRQ (ierr);

  /* Get pointer to vector data */
  ierr = VecGetArray(localX,&x); CHKERRQ (ierr);

  /* Get local grid boundaries */
  ierr = DAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm); CHKERRQ (ierr);
  ierr = DAGetGhostCorners(user->da,&gxs,&gys,&gzs,&gxm,&gym,&gzm);
  CHKERRQ (ierr);

  /* Determine the interior of the locally owned part of the grid. */
  if (user->period == DA_XYZPERIODIC) {
    xints = xs; xinte = xs+xm;
    yints = ys; yinte = ys+ym;
    zints = zs; zinte = zs+zm; }
  else {
    SETERRQ(PETSC_ERR_ARG_WRONGSTATE,
	    "Only periodic boundary conditions in 3-D Cahn-Hilliard"); }

  /* Loop over the interior points... */
  for (k=zints-gzs; k<zinte-gzs; k++)
    for (j=k*gym + yints-gys; j<k*gym + yinte-gys; j++)
      for (i=j*gxm + xints-gxs; i<j*gxm + xinte-gxs; i++) {

	/* Form the column indices */
	columns[0]=C(i+gxm*gym); columns[1]=C(i+gxm);
	columns[2]=C(i-1); columns[3]=C(i); columns[4]=C(i+1);
	columns[5]=C(i-gxm); columns[6]=C(i-gxm*gym);

	bvalue[0] = beta * hzm2 * ch_psidoubleprime (x[columns[0]], mparam);
	bvalue[1] = beta * hym2 * ch_psidoubleprime (x[columns[1]], mparam);
	bvalue[2] = beta * hxm2 * ch_psidoubleprime (x[columns[2]], mparam);
	bvalue[3] = -2.*beta * (hxm2+hym2+hzm2) *
	  ch_psidoubleprime (x[columns[3]], mparam);
	bvalue[4] = beta * hxm2 * ch_psidoubleprime (x[columns[4]], mparam);
	bvalue[5] = beta * hym2 * ch_psidoubleprime (x[columns[5]], mparam);
	bvalue[6] = beta * hzm2 * ch_psidoubleprime (x[columns[6]], mparam);

	node = C(i);
	MatSetValuesLocal (user->J, 1,&node, 5,columns, bvalue, ADD_VALUES);
      }

  /* Assemble matrix, using the 2-step process:
     MatAssemblyBegin(), MatAssemblyEnd(). */
  ierr = MatAssemblyBegin (user->J,MAT_FINAL_ASSEMBLY); CHKERRQ (ierr);
  ierr = MatAssemblyEnd (user->J,MAT_FINAL_ASSEMBLY); CHKERRQ (ierr);

  /* Restore unknown vector */
  ierr = VecRestoreArray (localX,&x); CHKERRQ (ierr);

  /* ierr = MatView (user->J, VIEWER_STDOUT_SELF); CHKERRQ (ierr); */

  return ierr;
}


/*++++++++++++++++++++++++++++++++++++++
  This creates the initial alpha and J matrices, where alpha is the alpha term
  component of the Jacobian.  Since the alpha term is linear, this part of the
  Jacobian need only be calculated once.

  int ch_jacobian_alpha_2d It returns zero or an error message.

  AppCtx *user The application context structure pointer.
  ++++++++++++++++++++++++++++++++++++++*/

int ch_jacobian_alpha_2d (AppCtx *user)
{
  int     ierr,node,i,j, mc,chvar, mx,my,xs,ys,xm,ym,gxs,gys,gxm,gym;
  int     xints,xinte,yints,yinte;
  int     columns [13];
  PetscScalar  hx,hy,dhx,dhy,hxm2,hym2;
  PetscScalar  alpha,beta,mparam;
  PetscScalar  avalue [25];

  mc = user->mc;  chvar = user->chvar;
  mx = user->mx;  my = user->my;  mparam = user->mparam;
  alpha = user->kappa * user->gamma * user->epsilon;
  beta = user->kappa * user->gamma / user->epsilon;

  /* Define mesh intervals ratios for uniform grid.
     [Note: FD formulae below may someday be normalized by multiplying through
     by local volume element to obtain coefficients O(1) in two dimensions.] */
  dhx = (PetscScalar)(mx-1);  dhy = (PetscScalar)(my-1);
  hx = 1./dhx;           hy = 1./dhy; /* This line hard-codes the 1x1 square */
  hxm2 = 1./hx/hx;       hym2 = 1./hy/hy;

  /* Get local grid boundaries */
  ierr = DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
  CHKERRQ (ierr);
  ierr = DAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL);
  CHKERRQ (ierr);

  /* Determine the interior of the locally owned part of the grid. */
  if (user->period == DA_XPERIODIC || user->period == DA_XYPERIODIC) {
    xints = xs; xinte = xs+xm; }
  else {
    xints = (xs==0) ? 2:xs; xinte = (xs+xm==mx) ? xs+xm-2 : xs+xm; }
  if (user->period == DA_YPERIODIC || user->period == DA_XYPERIODIC) {
    yints = ys; yinte = ys+ym; }
  else {
    yints = (ys==0) ? 2:ys; yinte = (ys+ym==my) ? ys+ym-2 : ys+ym; }

  /* Get parameters for matrix creation */
  i = xm * ym * user->mc;
  j = mx * my * user->mc;

  /* Create the distributed alpha matrix */
  ierr = MatCreateMPIAIJ (user->comm, i,i, j,j, 13,PETSC_NULL, 6,PETSC_NULL,
			  &(user->alpha));
  CHKERRQ (ierr);

  /* Get the local-to-global mapping and associate it with the matrix */
  ierr = DAGetISLocalToGlobalMapping (user->da, &user->isltog); CHKERRQ (ierr);
  ierr = MatSetLocalToGlobalMapping (user->alpha, user->isltog); CHKERRQ (ierr);

  /* Pre-compute alpha term values (see ch_residual_2d above)
     This should be the only place they're actually computed; they'll be
     used below. */
  avalue[0] = avalue[12] = -alpha*hym2*hym2;
  avalue[1] = avalue[3] = avalue[9] = avalue[11] = -2.*alpha*hxm2*hym2;
  avalue[2] = avalue[10] = 4.*alpha*hym2*(hxm2+hym2);
  avalue[4] = avalue[8] = -alpha*hxm2*hxm2;
  avalue[5] = avalue[7] = 4.*alpha*hxm2*(hxm2+hym2);
  avalue[6] = -alpha*(6.*hxm2*hxm2 + 6.*hym2*hym2 + 8.*hxm2*hym2);

  /* Loop over interior nodes */
  for (j=yints-gys; j<yinte-gys; j++)
    for (i=j*gxm + xints-gxs; i<j*gxm + xinte-gxs; i++) {

      /* Form the column indices */
      columns[0]=C(i+2*gxm);
      columns[1]=C(i+gxm-1); columns[2]=C(i+gxm);  columns[3]=C(i+gxm+1);
      columns[4]=C(i-2);     columns[5]=C(i-1);    columns[6]=C(i);
      columns[7]=C(i+1);     columns[8]=C(i+2);
      columns[9]=C(i-gxm-1); columns[10]=C(i-gxm); columns[11]=C(i-gxm+1);
      columns[12]=C(i-2*gxm);

      node = C(i);
      MatSetValuesLocal (user->alpha, 1,&node, 13,columns, avalue,
			 INSERT_VALUES);
    }

  /* Okay, that was the easy part, now for the boundary conditions! */

  /* Test whether we're on the bottom row and non-y-periodic */
  if (yints == 2) {
    printf ("This must only be reached if we are NOT y-periodic\n");
    /* Change value 6 and do the second-from-bottom row */
    avalue[6] += avalue[12];
    for (i=gxm + xints-gxs; i<gxm + xinte-gxs; i++) {

      /* Form the column indices */
      columns[0]=C(i+2*gxm);
      columns[1]=C(i+gxm-1); columns[2]=C(i+gxm); columns[3]=C(i+gxm+1);
      columns[4]=C(i-2);     columns[5]=C(i-1);   columns[6]=C(i);
      columns[7]=C(i+1);     columns[8]=C(i+2);
      columns[9]=C(i-gxm-1); columns[10]=C(i-gxm); columns[11]=C(i-gxm+1);

      node = C(i);
      MatSetValuesLocal (user->alpha, 1,&node, 12,columns, avalue,
			 INSERT_VALUES);
    }
    avalue[6] -= avalue[12];

    /* Make some new avalues and do the bottom row */
    avalue[13] = 2.*avalue[0];
    avalue[14] = avalue[16] = 2.*avalue[1];
    avalue[15] = 2.*avalue[2];
    avalue[17] = avalue[21] = avalue[4];
    avalue[18] = avalue[20] = avalue[5];
    avalue[19] = avalue[6];
    for (i=xints-gxs; i<xinte-gxs; i++) {

      /* Form the column indices */
      columns[0]=C(i+2*gxm);
      columns[1]=C(i+gxm-1); columns[2]=C(i+gxm); columns[3]=C(i+gxm+1);
      columns[4]=C(i-2);     columns[5]=C(i-1);   columns[6]=C(i);
      columns[7]=C(i+1);     columns[8]=C(i+2);

      node = C(i);
      MatSetValuesLocal (user->alpha, 1,&node, 9,columns, avalue+13,
			 INSERT_VALUES);
    }

    /* Bottom left corner */
    if (xints == 2) {
      node=C(0); /* The point (0,0) */
      columns[0]=C(2*gxm);
      columns[1]=C(gxm);   columns[2]=C(gxm+1);
      columns[3]=C(0);     columns[4]=C(1);     columns[5]=C(2);
      avalue[13] = 2.*avalue[0];
      avalue[14] = 2.*avalue[2];
      avalue[15] = 4.*avalue[3];
      avalue[16] = avalue[6];
      avalue[17] = 2.*avalue[7];
      avalue[18] = 2.*avalue[8];
      MatSetValuesLocal (user->alpha, 1,&node, 6,columns, avalue+13,
			 INSERT_VALUES);
      node=C(1); /* The point (1,0) */
      columns[0]=C(2*gxm+1);
      columns[1]=C(gxm); columns[2]=C(gxm+1); columns[3]=C(gxm+2);
      columns[4]=C(0);   columns[5]=C(1);  columns[6]=C(2);  columns[7]=C(3);
      avalue[13] = 2.*avalue[0];
      avalue[14] = avalue[16] = 2.*avalue[1];
      avalue[15] = 2.*avalue[2];
      avalue[17] = avalue[19] = avalue[5];
      avalue[18] = avalue[6]  + avalue[4];
      avalue[20] = avalue[8];
      MatSetValuesLocal (user->alpha, 1,&node, 8,columns, avalue+13,
			 INSERT_VALUES);
      node=C(gxm); /* The point (0,1) */
      columns[0]=C(3*gxm);
      columns[1]=C(2*gxm); columns[2]=C(2*gxm+1);
      columns[3]=C(gxm);   columns[4]=C(gxm+1);   columns[5]=C(gxm+2);
      columns[6]=C(0);     columns[7]=C(1);
      avalue[13] = avalue[0];
      avalue[14] = avalue[19] = avalue[2];
      avalue[15] = avalue[20] = 2.*avalue[3];
      avalue[16] = avalue[6] + avalue[12];
      avalue[17] = 2.*avalue[7];
      avalue[18] = 2.*avalue[8];
      MatSetValuesLocal (user->alpha, 1,&node, 8,columns, avalue+13,
			 INSERT_VALUES);
      node=C(gxm+1); /* The point (1,1) */
      columns[0]=C(3*gxm+1);
      columns[1]=C(2*gxm); columns[2]=C(2*gxm+1); columns[3]=C(2*gxm+2);
      columns[4]=C(gxm);   columns[5]=C(gxm+1);   columns[6]=C(gxm+2);
      columns[7]=C(gxm+3);
      columns[8]=C(0);     columns[9]=C(1);       columns[10]=C(2);
      avalue[13] = avalue[0];
      avalue[14] = avalue[16] = avalue[21] = avalue[23] = avalue[1];
      avalue[15] = avalue[22] = avalue[2];
      avalue[17] = avalue[19] = avalue[5];
      avalue[18] = avalue[6] + avalue[4] + avalue[12];
      avalue[20] = avalue[8];
      MatSetValuesLocal (user->alpha, 1,&node, 11,columns, avalue+13,
			 INSERT_VALUES);
    }

    /* Bottom right corner */
    if(xinte == mx-2) {
      node=C(i=mx-gxs-1); /* The point (mx-1, 0) */
      columns[0]=C(i+2*gxm);
      columns[1]=C(i+gxm-1); columns[2]=C(i+gxm);
      columns[3]=C(i-2);     columns[4]=C(i-1);   columns[5]=C(i);
      avalue[13] = 2.*avalue[0];
      avalue[14] = 4.*avalue[1];
      avalue[15] = 2.*avalue[2];
      avalue[16] = 2.*avalue[4];
      avalue[17] = 2.*avalue[5];
      avalue[18] = avalue[6];
      MatSetValuesLocal (user->alpha, 1,&node, 6,columns, avalue+13,
			 INSERT_VALUES);
      node=C(i=mx-gxs-2); /* The point (mx-2,0) */
      columns[0]=C(i+2*gxm);
      columns[1]=C(i+gxm-1); columns[2]=C(i+gxm); columns[3]=C(i+gxm+1);
      columns[4]=C(i-2);     columns[5]=C(i-1);   columns[6]=C(i);
      columns[7]=C(i+1);
      avalue[13] = 2.*avalue[0];
      avalue[14] = avalue[16] = 2.*avalue[1];
      avalue[15] = 2.*avalue[2];
      avalue[17] = avalue[4];
      avalue[18] = avalue[20] = avalue[5];
      avalue[19] = avalue[6]  + avalue[8];
      MatSetValuesLocal (user->alpha, 1,&node, 8,columns, avalue+13,
			 INSERT_VALUES);
      node=C(i=gxm+mx-gxs-1); /* The point (mx-1,1) */
      columns[0]=C(i+2*gxm);
      columns[1]=C(i+gxm-1); columns[2]=C(i+gxm);
      columns[3]=C(i-2);     columns[4]=C(i-1);   columns[5]=C(i);
      columns[6]=C(i-gxm-1); columns[7]=C(i-gxm);
      avalue[13] = avalue[0];
      avalue[14] = avalue[19] = 2.*avalue[1];
      avalue[15] = avalue[20] = avalue[2];
      avalue[16] = 2.*avalue[4];
      avalue[17] = 2.*avalue[5];
      avalue[18] = avalue[6] + avalue[12];
      MatSetValuesLocal (user->alpha, 1,&node, 8,columns, avalue+13,
			 INSERT_VALUES);
      node=C(i=gxm+mx-gxs-2); /* The point (mx-2,1) */
      columns[0]=C(i+2*gxm);
      columns[1]=C(i+gxm-1); columns[2]=C(i+gxm); columns[3]=C(i+gxm+1);
      columns[4]=C(i-2);     columns[5]=C(i-1);   columns[6]=C(i);
      columns[7]=C(i+1);
      columns[8]=C(i-gxm-1); columns[9]=C(i-gxm);   columns[10]=C(i-gxm+1);
      avalue[13] = avalue[0];
      avalue[14] = avalue[16] = avalue[21] = avalue[23] = avalue[1];
      avalue[15] = avalue[22] = avalue[2];
      avalue[17] = avalue[4];
      avalue[18] = avalue[20] = avalue[5];
      avalue[19] = avalue[6] + avalue[8] + avalue[12];
      MatSetValuesLocal (user->alpha, 1,&node, 11,columns, avalue+13,
			 INSERT_VALUES);
    }
  }

  /* Test whether we're on the top row and non-y-periodic */
  if (yinte == my-2) {
    printf ("This must only be reached if we are NOT y-periodic\n");
    /* Change value 6 and do the second-from-top row */
    avalue[6] += avalue[0];
    for (i=(my-gys-2)*gxm + xints-gxs;
	 i<(my-gys-2)*gxm + xinte-gxs; i++) {

      /* Form the column indices */
      columns[1]=C(i+gxm-1); columns[2]=C(i+gxm);  columns[3]=C(i+gxm+1);
      columns[4]=C(i-2);     columns[5]=C(i-1);    columns[6]=C(i);
      columns[7]=C(i+1);     columns[8]=C(i+2);
      columns[9]=C(i-gxm-1); columns[10]=C(i-gxm); columns[11]=C(i-gxm+1);
      columns[12]=C(i-2*gxm);

      node = C(i);
      MatSetValuesLocal (user->alpha, 1,&node, 12,columns+1, avalue+1,
			 INSERT_VALUES);
    }
    avalue[6] -= avalue[0];

    /* Make some new avalues and do the top row */
    avalue[13] = avalue[17] = avalue[4];
    avalue[14] = avalue[16] = avalue[5];
    avalue[15] = avalue[6];
    avalue[18] = avalue[20] = 2.*avalue[9];
    avalue[19] = 2.*avalue[10];
    avalue[21] = 2.*avalue[12];
    for (i=(my-gys-1)*gxm + xints-gxs;
	 i<(my-gys-1)*gxm + xinte-gxs; i++) {

      /* Form the column indices */
      columns[0]=C(i-2);     columns[1]=C(i-1);   columns[2]=C(i);
      columns[3]=C(i+1);     columns[4]=C(i+2);
      columns[5]=C(i-gxm-1); columns[6]=C(i-gxm); columns[7]=C(i-gxm+1);
      columns[8]=C(i-2*gxm);

      node = C(i);
      MatSetValuesLocal (user->alpha, 1,&node, 9,columns, avalue+13,
			 INSERT_VALUES);
    }

    /* Top left corner */
    if (xints == 2) {
      node=C(i=(my-gys-1)*gxm); /* The point (0,my-1) */
      columns[0]=C(i);       columns[1]=C(i+1); columns[2]=C(i+2);
      columns[3]=C(i-gxm);   columns[4]=C(i-gxm+1);
      columns[5]=C(i-2*gxm);
      avalue[13] = avalue[6];
      avalue[14] = 2.*avalue[7];
      avalue[15] = 2.*avalue[8];
      avalue[16] = 2.*avalue[10];
      avalue[17] = 4.*avalue[11];
      avalue[18] = 2.*avalue[12];
      MatSetValuesLocal (user->alpha, 1,&node, 6,columns, avalue+13,
			 INSERT_VALUES);
      node=C(i=(my-gys-1)*gxm+1); /* The point (1,my-1) */
      columns[0]=C(i-1);     columns[1]=C(i);     columns[2]=C(i+1);
      columns[3]=C(i+2);
      columns[4]=C(i-gxm-1); columns[5]=C(i-gxm); columns[6]=C(i-gxm+1);
      columns[7]=C(i-2*gxm);
      avalue[13] = avalue[15] = avalue[5];
      avalue[14] = avalue[6]  + avalue[4];
      avalue[16] = avalue[8];
      avalue[17] = avalue[19] = 2.*avalue[9];
      avalue[18] = 2.*avalue[10];
      avalue[20] = 2.*avalue[12];
      MatSetValuesLocal (user->alpha, 1,&node, 8,columns, avalue+13,
			 INSERT_VALUES);
      node=C(i=(my-gys-2)*gxm); /* The point (0,my-2) */
      columns[0]=C(i+gxm); columns[1]=C(i+gxm+1);
      columns[2]=C(i);     columns[3]=C(i+1);   columns[4]=C(i+2);
      columns[5]=C(i-gxm); columns[6]=C(i-gxm+1);
      columns[7]=C(i-2*gxm);
      avalue[13] = avalue[18] = avalue[2];
      avalue[14] = avalue[19] = 2.*avalue[3];
      avalue[15] = avalue[6] + avalue[0];
      avalue[16] = 2.*avalue[7];
      avalue[17] = 2.*avalue[8];
      avalue[20] = avalue[12];
      MatSetValuesLocal (user->alpha, 1,&node, 8,columns, avalue+13,
			 INSERT_VALUES);
      node=C(i=(my-gys-2)*gxm+1); /* The point (1,my-2) */
      columns[0]=C(i+gxm-1); columns[1]=C(i+gxm); columns[2]=C(i+gxm+1);
      columns[3]=C(i-1);     columns[4]=C(i);     columns[5]=C(i+1);
      columns[6]=C(i+2);
      columns[7]=C(i-gxm-1); columns[8]=C(i-gxm); columns[9]=C(i-gxm+1);
      columns[10]=C(i-2*gxm);
      avalue[13] = avalue[15] = avalue[20] = avalue[22] = avalue[1];
      avalue[14] = avalue[21] = avalue[2];
      avalue[16] = avalue[18] = avalue[5];
      avalue[17] = avalue[6] + avalue[4] + avalue[0];
      avalue[19] = avalue[8];
      avalue[23] = avalue[12];
      MatSetValuesLocal (user->alpha, 1,&node, 11,columns, avalue+13,
			 INSERT_VALUES);
    }

    /* Top right corner */
    if(xinte == mx-2) {
      node=C(i=(my-gys-1)*gxm + mx-gxs-1); /* The point (mx-1, my-1) */
      columns[0]=C(i-2);   columns[1]=C(i-1);   columns[2]=C(i);
      columns[3]=C(i-gxm-1); columns[4]=C(i-gxm);
      columns[5]=C(i-2*gxm);
      avalue[13] = 2.*avalue[4];
      avalue[14] = 2.*avalue[5];
      avalue[15] = avalue[6];
      avalue[16] = 4.*avalue[9];
      avalue[17] = 2.*avalue[10];
      avalue[18] = 2.*avalue[12];
      MatSetValuesLocal (user->alpha, 1,&node, 6,columns, avalue+13,
			 INSERT_VALUES);
      node=C(i=(my-gys-1)*gxm + mx-gxs-2); /* The point (mx-2, my-1) */
      columns[0]=C(i-2);   columns[1]=C(i-1);   columns[2]=C(i);
      columns[3]=C(i+1);
      columns[4]=C(i-gxm-1); columns[5]=C(i-gxm); columns[6]=C(i-gxm+1);
      columns[7]=C(i-2*gxm);
      avalue[13] = avalue[4];
      avalue[14] = avalue[16] = avalue[5];
      avalue[15] = avalue[6]  + avalue[8];
      avalue[17] = avalue[19] = 2.*avalue[9];
      avalue[18] = 2.*avalue[10];
      avalue[20] = 2.*avalue[12];
      MatSetValuesLocal (user->alpha, 1,&node, 8,columns, avalue+13,
			 INSERT_VALUES);
      node=C(i=(my-gys-2)*gxm + mx-gxs-1); /* The point (mx-1, my-2) */
      columns[0]=C(i+gxm-1); columns[1]=C(i+gxm);
      columns[2]=C(i-2);     columns[3]=C(i-1);   columns[4]=C(i);
      columns[5]=C(i-gxm-1); columns[6]=C(i-gxm);
      columns[7]=C(i-2*gxm);
      avalue[13] = avalue[18] = 2.*avalue[1];
      avalue[14] = avalue[19] = avalue[2];
      avalue[15] = 2.*avalue[4];
      avalue[16] = 2.*avalue[5];
      avalue[17] = avalue[6] + avalue[0];
      avalue[20] = avalue[12];
      MatSetValuesLocal (user->alpha, 1,&node, 8,columns, avalue+13,
			 INSERT_VALUES);
      node=C(i=(my-gys-2)*gxm + mx-gxs-2); /* The point (mx-2, my-2) */
      columns[1]=C(i+gxm-1); columns[2]=C(i+gxm); columns[3]=C(i+gxm+1);
      columns[4]=C(i-2);     columns[5]=C(i-1);   columns[6]=C(i);
      columns[7]=C(i+1);
      columns[8]=C(i-gxm-1); columns[9]=C(i-gxm); columns[10]=C(i-gxm+1);
      columns[11]=C(i-2*gxm);
      avalue[14] = avalue[16] = avalue[21] = avalue[23] = avalue[1];
      avalue[15] = avalue[22] = avalue[2];
      avalue[17] = avalue[4];
      avalue[18] = avalue[20] = avalue[5];
      avalue[19] = avalue[6] + avalue[8] + avalue[0];
      avalue[24] = avalue[12];
      MatSetValuesLocal (user->alpha, 1,&node, 11,columns+1, avalue+14,
			 INSERT_VALUES);
    }
  }

  /* Test whether we're on the left column and non-x-periodic */
  if (xints == 2) {
    printf ("This must only be reached if we are NOT x-periodic\n");
    /* Make some avalues and do the second-from-left column */
    avalue[13] = avalue[24] = avalue[0];
    avalue[14] = avalue[16] = avalue[21] = avalue[23] = avalue[1];
    avalue[15] = avalue[22] = avalue[2];
    avalue[17] = avalue[19] = avalue[5];
    avalue[18] = avalue[6]  + avalue[4];
    avalue[20] = avalue[8];
    for (i=(yints-gys)*gxm + 1; i<(yinte-gys)*gxm + 1; i+=gxm) {
      columns[0]=C(i+2*gxm);
      columns[1]=C(i+gxm-1); columns[2]=C(i+gxm);  columns[3]=C(i+gxm+1);
      columns[4]=C(i-1);     columns[5]=C(i);      columns[6]=C(i+1);
      columns[7]=C(i+2);
      columns[8]=C(i-gxm-1); columns[9]=C(i-gxm); columns[10]=C(i-gxm+1);
      columns[11]=C(i-2*gxm);

      node = C(i);
      MatSetValuesLocal (user->alpha, 1,&node, 12,columns, avalue+13,
			 INSERT_VALUES);
    }

    /* Make some more avalues and do the leftmost column */
    avalue[13] = avalue[21] = avalue[0];
    avalue[14] = avalue[19] = avalue[2];
    avalue[15] = avalue[20] = 2.*avalue[3];
    avalue[16] = avalue[6];
    avalue[17] = 2.*avalue[7];
    avalue[18] = 2.*avalue[8];
    for (i=(yints-gys)*gxm; i<(yinte-gys)*gxm; i+=gxm) {
      columns[0]=C(i+2*gxm);
      columns[1]=C(i+gxm);  columns[2]=C(i+gxm+1);
      columns[3]=C(i);      columns[4]=C(i+1);     columns[5]=C(i+2);
      columns[6]=C(i-gxm);  columns[7]=C(i-gxm+1);
      columns[8]=C(i-2*gxm);

      node = C(i);
      MatSetValuesLocal (user->alpha, 1,&node, 9,columns, avalue+13,
			 INSERT_VALUES);
    }
  }

  /* Test whether we're on the right column and non-x-periodic */
  if (xinte == mx-2) {
    printf ("This must only be reached if we are NOT x-periodic\n");
    /* Make some avalues and do the second-from-right column */
    avalue[13] = avalue[24] = avalue[0];
    avalue[14] = avalue[16] = avalue[21] = avalue[23] = avalue[1];
    avalue[15] = avalue[22] = avalue[2];
    avalue[17] = avalue[4];
    avalue[18] = avalue[20] = avalue[5];
    avalue[19] = avalue[6]  + avalue[8];
    for (i=(yints-gys)*gxm + mx-gxs-2; i<(yinte-gys)*gxm + mx-gxs-2; i+=gxm) {
      columns[0]=C(i+2*gxm);
      columns[1]=C(i+gxm-1); columns[2]=C(i+gxm);  columns[3]=C(i+gxm+1);
      columns[4]=C(i-2);     columns[5]=C(i-1);    columns[6]=C(i);
      columns[7]=C(i+1);
      columns[8]=C(i-gxm-1); columns[9]=C(i-gxm); columns[10]=C(i-gxm+1);
      columns[11]=C(i-2*gxm);

      node = C(i);
      MatSetValuesLocal (user->alpha, 1,&node, 12,columns, avalue+13,
			 INSERT_VALUES);
    }

    /* Make some more avalues and do the rightmost column */
    avalue[13] = avalue[21] = avalue[0];
    avalue[14] = avalue[19] = 2.*avalue[1];
    avalue[15] = avalue[20] = avalue[2];
    avalue[16] = 2.*avalue[4];
    avalue[17] = 2.*avalue[5];
    avalue[18] = avalue[6];
    for (i=(yints-gys)*gxm + mx-gxs-1; i<(yinte-gys)*gxm + mx-gxs-1; i+=gxm) {
      columns[0]=C(i+2*gxm);
      columns[1]=C(i+gxm-1);  columns[2]=C(i+gxm);
      columns[3]=C(i-2);      columns[4]=C(i-1);     columns[5]=C(i);
      columns[6]=C(i-gxm-1);  columns[7]=C(i-gxm);
      columns[8]=C(i-2*gxm);

      node = C(i);
      MatSetValuesLocal (user->alpha, 1,&node, 9,columns, avalue+13,
			 INSERT_VALUES);
    }
  }

  /* Assemble matrix, using the 2-step process:
     MatAssemblyBegin(), MatAssemblyEnd(). */
  ierr = MatAssemblyBegin(user->alpha,MAT_FINAL_ASSEMBLY); CHKERRQ (ierr);
  ierr = MatAssemblyEnd(user->alpha,MAT_FINAL_ASSEMBLY); CHKERRQ (ierr);

  /* Make J a copy of alpha with the same local mapping (setting new mapping is
     unnecessary with PETSc 2.1.5). */
  ierr = MatDuplicate (user->alpha, MAT_COPY_VALUES, &(user->J)); CHKERRQ (ierr);

  /* Tell the matrix we will never add a new nonzero location to the
     matrix. If we do, it will generate an error. */
  ierr = MatSetOption(user->J,MAT_NEW_NONZERO_LOCATION_ERR); CHKERRQ (ierr);

  return ierr;
}


/*++++++++++++++++++++++++++++++++++++++
  This creates the initial alpha and J matrices, where alpha is the alpha term
  component of the Jacobian.  Since the alpha term is linear, this part of the
  Jacobian need only be calculated once.

  int ch_jacobian_alpha_3d It returns zero or an error message.

  AppCtx *user The application context structure pointer.
  ++++++++++++++++++++++++++++++++++++++*/

int ch_jacobian_alpha_3d (AppCtx *user)
{
  int     ierr,node,i,j,k, mc,chvar, mx,my,mz, xs,ys,zs, xm,ym,zm;
  int     gxs,gys,gzs, gxm,gym,gzm, xints,xinte, yints,yinte, zints,zinte;
  int     columns [25];
  PetscScalar  hx,hy,hz, dhx,dhy,dhz, hxm2,hym2,hzm2;
  PetscScalar  alpha,beta,mparam;
  PetscScalar  avalue [25];

  mc = user->mc;  chvar = user->chvar;
  mx = user->mx;  my = user->my;  mz = user->mz;  mparam = user->mparam;
  alpha = user->kappa * user->gamma * user->epsilon;
  beta = user->kappa * user->gamma / user->epsilon;

  /* Define mesh intervals ratios for uniform grid.
     [Note: FD formulae below may someday be normalized by multiplying through
     by local volume element to obtain coefficients O(1) in two dimensions.] */
  dhx = (PetscScalar)(mx-1);  dhy = (PetscScalar)(my-1);
  dhz = (PetscScalar)(mz-1);
  /* This line hard-codes the 1x1x1 cube */
  hx = 1./dhx;           hy = 1./dhy;          hz = 1./dhz;
  hxm2 = 1./hx/hx;       hym2 = 1./hy/hy;      hzm2 = 1./hz/hz;

  /* Get local grid boundaries */
  ierr = DAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm); CHKERRQ (ierr);
  ierr = DAGetGhostCorners(user->da,&gxs,&gys,&gzs,&gxm,&gym,&gzm);
  CHKERRQ (ierr);

  /* Determine the interior of the locally owned part of the grid. */
  if (user->period == DA_XYZPERIODIC) {
    xints = xs; xinte = xs+xm;
    yints = ys; yinte = ys+ym;
    zints = zs; zinte = zs+zm; }
  else {
    SETERRQ(PETSC_ERR_ARG_WRONGSTATE,
	    "Only periodic boundary conditions in 3-D Cahn-Hilliard"); }

  /* Get parameters for matrix creation */
  i = xm * ym * zm * user->mc;
  j = mx * my * mz * user->mc;

  /* Create the distributed alpha matrix */
  ierr = MatCreateMPIAIJ (user->comm, i,i, j,j, 25,PETSC_NULL, 12,PETSC_NULL,
			  &(user->alpha));
  CHKERRQ (ierr);

  /* Get the local-to-global mapping and associate it with the matrix */
  ierr = DAGetISLocalToGlobalMapping (user->da, &user->isltog); CHKERRQ (ierr);
  ierr = MatSetLocalToGlobalMapping (user->alpha, user->isltog); CHKERRQ (ierr);

  /* Pre-compute alpha term values (see ch_residual_2d above)
     This should be the only place they're actually computed; they'll be
     used below. */
  avalue[0]  = avalue[24] = -alpha*hzm2*hzm2;
  avalue[1]  = avalue[5] = avalue[19] = avalue[23] = -2.*alpha*hym2*hzm2;
  avalue[2]  = avalue[4] = avalue[20] = avalue[22] = -2.*alpha*hxm2*hzm2;
  avalue[3]  = avalue[21] = 4.*alpha*hzm2*(hxm2+hym2+hzm2);
  avalue[6]  = avalue[18] = -alpha*hym2*hym2;
  avalue[7]  = avalue[9] = avalue[15] = avalue[17] = -2.*alpha*hxm2*hym2;
  avalue[8]  = avalue[16] = 4.*alpha*hym2*(hxm2+hym2+hzm2);
  avalue[10] = avalue[14] = -alpha*hxm2*hxm2;
  avalue[11] = avalue[13] = 4.*alpha*hxm2*(hxm2+hym2+hzm2);
  avalue[12] = -alpha*(6.*hxm2*hxm2 + 6.*hym2*hym2 + 6.*hzm2*hzm2
		       + 8.*hxm2*hym2 + 8.*hxm2*hzm2 + 8.*hym2*hzm2);

  /* Loop over interior nodes */
  for (k=zints-gzs; k<zinte-gzs; k++)
    for (j=k*gym + yints-gys; j<k*gym + yinte-gys; j++)
      for (i=j*gxm + xints-gxs; i<j*gxm + xinte-gxs; i++) {

	/* Form the column indices */
	columns[0]=C(i+2*gxm*gym);
	columns[1]=C(i+gxm*gym+gxm);
	columns[2]=C(i+gxm*gym-1);
	columns[3]=C(i+gxm*gym);
	columns[4]=C(i+gxm*gym+1);
	columns[5]=C(i+gxm*gym-gxm);
	columns[6]=C(i+2*gxm);
	columns[7]=C(i+gxm-1);  columns[8]=C(i+gxm);  columns[9]=C(i+gxm+1);
	columns[10]=C(i-2);     columns[11]=C(i-1);   columns[12]=C(i);
	columns[13]=C(i+1);     columns[14]=C(i+2);
	columns[15]=C(i-gxm-1); columns[16]=C(i-gxm); columns[17]=C(i-gxm+1);
	columns[18]=C(i-2*gxm);
	columns[19]=C(i-gxm*gym+gxm);
	columns[20]=C(i-gxm*gym-1);
	columns[21]=C(i-gxm*gym);
	columns[22]=C(i-gxm*gym+1);
	columns[23]=C(i-gxm*gym-gxm);
	columns[24]=C(i-2*gxm*gym);

	node = C(i);
	MatSetValuesLocal (user->alpha, 1,&node, 25,columns, avalue,
			   INSERT_VALUES);
      }

  /* Assemble matrix, using the 2-step process:
     MatAssemblyBegin(), MatAssemblyEnd(). */
  ierr = MatAssemblyBegin(user->alpha,MAT_FINAL_ASSEMBLY); CHKERRQ (ierr);
  ierr = MatAssemblyEnd(user->alpha,MAT_FINAL_ASSEMBLY); CHKERRQ (ierr);

  /* Make J a copy of alpha with the same local mapping (setting new mapping is
     unnecessary with PETSc 2.1.5). */
  ierr = MatDuplicate (user->alpha, MAT_COPY_VALUES, &(user->J)); CHKERRQ (ierr);

  /* Tell the matrix we will never add a new nonzero location to the
     matrix. If we do, it will generate an error. */
  ierr = MatSetOption(user->J,MAT_NEW_NONZERO_LOCATION_ERR); CHKERRQ (ierr);

  return ierr;
}