1    | /***************************************
2    |   $Header$
3    | 
4    |   This file contains the heart of the Cahn-Hilliard formulation, in particular
5    |   the functions which build the finite difference residuals and Jacobian.
6    | ***************************************/
7    | 
8    | 
9    | /* Including cahnhill.h includes the necessary PETSc header files */
10   | #include "cahnhill.h"
11   | 
12   | 
13   | /*++++++++++++++++++++++++++++++++++++++
14   |   This calculates the derivative of homogeneous free energy with respect to
15   |   concentration, which is a component of the chemical potential.  It currently
16   |   uses
17   |   +latex+$$\Psi' = C(1-C)\left(\frac{1}{2}+m-C\right)$$
18   |   +html+ <center>Psi' = <i>C</i> (1-<i>C</i>) (0.5+<i>m</i>+<i>C</i>)</center>
19   |   which gives (meta)stable equilibria at
20   |   +latex+$C=0$ and 1 and an unstable equilibrium at $C=\frac{1}{2}+m$; if $m>0$
21   |   +html+ <i>C</i>=0 and 1 and an unstable equilibrium at <i>C</i> = 1/2 +
22   |   +html+ <i>m</i>; if <i>m</i>&gt;0
23   |   then the 0 phase is stable and vice versa.
24   | 
25   |   PetscScalar ch_psiprime It returns the derivative itself.
26   | 
27   |   PetscScalar C The concentration.
28   | 
29   |   PetscScalar mparam The model parameter
30   |   +latex+$m$.
31   |   +html+ <i>m</i>.
32   |   ++++++++++++++++++++++++++++++++++++++*/
33   | 
34   | static inline PetscScalar ch_psiprime (PetscScalar C, PetscScalar mparam)
35   | { return C*(1-C)*(0.5+mparam-C); }
36   | 
37   | #define PSIPRIME_FLOPS 5
38   | 
39   | 
40   | /*++++++++++++++++++++++++++++++++++++++
41   |   This calculates the second derivative of homogeneous free energy with respect
42   |   to concentration, for insertion into the Jacobian matrix.  See the
43   |   +latex+{\tt ch\_psiprime()} function for the definition of $\Psi$.
44   |   +html+ <tt>ch_psiprime()</tt> function for the definition of Psi.
45   | 
46   |   PetscScalar ch_psidoubleprime It returns the second derivative
47   |   +latex+$\Psi''(C)$.
48   |   -latex-Psi''(C).
49   | 
50   |   PetscScalar C The concentration.
51   | 
52   |   PetscScalar mparam The model parameter
53   |   +latex+$m$.
54   |   +html+ <i>m</i>.
55   |   ++++++++++++++++++++++++++++++++++++++*/
56   | 
57   | static inline PetscScalar ch_psidoubleprime (PetscScalar C, PetscScalar mparam)
58   | { return 3*C*C - (3+2*mparam)*C + 0.5+mparam; }
59   | 
60   | #define PSIDOUBLEPRIME_FLOPS 8
61   | 
62   | 
63   | /*++++++++++++++++++++++++++++++++++++++
64   |   This function computes the residual from indices to points in the
65   |   concentration array.  ``Up'' refers to the positive
66   |   +latex+$y$-direction, ``down'' to negative $y$, ``left'' to negative $x$ and
67   |   +latex+``right'' to positive $x$.
68   |   +html+ <i>y</i>-direction, ``down'' to negative <i>y</i>, ``left'' to
69   |   +html+ negative <i>x</i> and ``right'' to positive <i>x</i>.
70   | 
71   |   PetscScalar ch_residual_2d Returns the residual itself
72   | 
73   |   PetscScalar *conc Array of concentrations
74   | 
75   |   PetscScalar alpha Model parameter
76   |   +latex+$\kappa\alpha$
77   |   -latex-kappa alpha
78   | 
79   |   PetscScalar beta Model parameter
80   |   +latex+$\kappa\beta$
81   |   -latex-kappa beta
82   | 
83   |   PetscScalar mparam Model parameter
84   |   +latex+$m$
85   |   -latex-m
86   | 
87   |   PetscScalar hx Inverse square
88   |   +latex+$x$-spacing $h_x^{-2}$
89   |   +html+ <i>x</i>-spacing <i>h<sub>x</sub></i><sup>-2</sup>
90   | 
91   |   PetscScalar hy Inverse square
92   |   +latex+$y$-spacing $h_y^{-2}$
93   |   +html+ <i>y</i>-spacing <i>h<sub>y</sub></i><sup>-2</sup>
94   | 
95   |   int upup Index to array position two cells up from current
96   | 
97   |   int upleft Index to array position one cell up and one left from current
98   | 
99   |   int up Index to array position one cell up from current
100  | 
101  |   int upright Index to array position one cell up and one right from current
102  | 
103  |   int leftleft Index to array position two cells left of current
104  | 
105  |   int left Index to array position one cell left of current
106  | 
107  |   int current Index to current cell array position
108  | 
109  |   int right Index to array position one cell right of current
110  | 
111  |   int rightright Index to array position two cells right of current
112  | 
113  |   int downleft Index to array position one cell down and one left from current
114  | 
115  |   int down Index to array position one cell down from current
116  | 
117  |   int downright Index to array position one cell down and one right from current
118  | 
119  |   int downdown Index to array position two cells down from current
120  |   ++++++++++++++++++++++++++++++++++++++*/
121  | 
122  | static inline PetscScalar ch_residual_2d
123  | (PetscScalar *conc, PetscScalar alpha, PetscScalar beta, PetscScalar mparam,
124  |  PetscScalar hx, PetscScalar hy, int upup, int upleft, int up, int upright, int
125  |  leftleft, int left, int current, int right, int rightright, int downleft, int
126  |  down, int downright, int downdown)
127  | {
128  |   PetscScalar retval, hx2, hy2, hxhy;
129  | 
130  |   hx2 = hx*hx;  hy2 = hy*hy;  hxhy = hx*hy;
131  | 
132  |   /*+ This calculates the
133  |     +latex+$\beta$-term, $\kappa\beta$ times the Laplacian of $\Psi'(C)$,
134  |     -latex-beta term, kappa*beta times the laplacian of Psi prime(C),
135  |     +*/
136  |   retval = beta *
137  |     (hx * (ch_psiprime(conc[left], mparam) + ch_psiprime(conc[right], mparam))
138  |      + hy * (ch_psiprime(conc[up], mparam) + ch_psiprime(conc[down], mparam))
139  |      - 2*(hx+hy) * ch_psiprime(conc[current], mparam));
140  | 
141  |   /*+ then subtracts the
142  |     +latex+$\alpha$-term, $\kappa\alpha\nabla^2\nabla^2C$.
143  |     -latex-alpha term, kappa*alpha times the Laplacian of the Laplacian of C.
144  |     +*/
145  |   retval -= alpha *
146  |     (hx2 * (conc[leftleft] + conc[rightright])
147  |      + hy2 * (conc[upup] + conc[downdown])
148  |      + 2*hxhy * (conc[upleft]+conc[upright]+conc[downleft]+conc[downright])
149  |      - (4*hx2 + 4*hxhy) * (conc[left] + conc[right])
150  |      - (4*hy2 + 4*hxhy) * (conc[up] + conc[down])
151  |      + (6*hx2 + 6*hy2 + 8*hxhy) * conc[current]);
152  | 
153  |   return retval;
154  | }
155  | 
156  | #define RESIDUAL_FLOPS_2D (5*PSIPRIME_FLOPS + 10 + 32)
157  | 
158  | 
159  | /*++++++++++++++++++++++++++++++++++++++
160  |   This function computes the residual from indices to points in the
161  |   concentration array.  ``Front refers to the positive
162  |   +latex+$z$-direction, ``back'' to negative $z$, ``up'' to positive
163  |   +latex+$y$, ``down'' to negative $y$, ``left'' to
164  |   +latex+negative $x$ and ``right'' to positive $x$.
165  |   +html+ <i>z</i>-direction, ``back'' to negative <i>z</i>, ``up'' to positive
166  |   +html+ <i>y</i>, ``down'' to negative <i>y</i>, ``left'' to
167  |   +html+ negative <i>x</i> and ``right'' to positive <i>x</i>.
168  | 
169  |   PetscScalar ch_residual_3d Returns the residual itself
170  | 
171  |   PetscScalar *conc Array of concentrations
172  | 
173  |   PetscScalar alpha Model parameter
174  |   +latex+$\kappa\alpha$
175  |   -latex-kappa alpha
176  | 
177  |   PetscScalar beta Model parameter
178  |   +latex+$\kappa\beta$
179  |   -latex-kappa beta
180  | 
181  |   PetscScalar mparam Model parameter
182  |   +latex+$m$
183  |   -latex-m
184  | 
185  |   PetscScalar hx Inverse square
186  |   +latex+$x$-spacing $h_x^{-2}$
187  |   +html+ <i>x</i>-spacing <i>h<sub>x</sub></i><sup>-2</sup>
188  | 
189  |   PetscScalar hy Inverse square
190  |   +latex+$y$-spacing $h_y^{-2}$
191  |   +html+ <i>y</i>-spacing <i>h<sub>y</sub></i><sup>-2</sup>
192  | 
193  |   PetscScalar hz Inverse square
194  |   +latex+$z$-spacing $h_z^{-2}$
195  |   +html+ <i>z</i>-spacing <i>h<sub>z</sub></i><sup>-2</sup>
196  | 
197  |   int frontfront Index to array position two cells in front of current
198  | 
199  |   int frontup Index to array position one cell front and one up from current
200  | 
201  |   int frontleft Index to array position one cell front and one left from current
202  | 
203  |   int front Index to array position one cell in front of current
204  | 
205  |   int frontright Index to array position one cell front and one right from current
206  | 
207  |   int frontdown Index to array position one cell front and one down from current
208  | 
209  |   int upup Index to array position two cells up from current
210  | 
211  |   int upleft Index to array position one cell up and one left from current
212  | 
213  |   int up Index to array position one cell up from current
214  | 
215  |   int upright Index to array position one cell up and one right from current
216  | 
217  |   int leftleft Index to array position two cells left of current
218  | 
219  |   int left Index to array position one cell left of current
220  | 
221  |   int current Index to current cell array position
222  | 
223  |   int right Index to array position one cell right of current
224  | 
225  |   int rightright Index to array position two cells right of current
226  | 
227  |   int downleft Index to array position one cell down and one left from current
228  | 
229  |   int down Index to array position one cell down from current
230  | 
231  |   int downright Index to array position one cell down and one right from current
232  | 
233  |   int downdown Index to array position two cells down from current
234  | 
235  |   int backup Index to array position one cell back and one up from current
236  | 
237  |   int backleft Index to array position one cell back and one left from current
238  | 
239  |   int back Index to array position one cell in back of current
240  | 
241  |   int backright Index to array position one cell back and one right from current
242  | 
243  |   int backdown Index to array position one cell back and one down from current
244  | 
245  |   int backback Index to array position two cells in back of current
246  |   ++++++++++++++++++++++++++++++++++++++*/
247  | 
248  | static inline PetscScalar ch_residual_3d
249  | (PetscScalar *conc, PetscScalar alpha, PetscScalar beta, PetscScalar mparam,
250  |  PetscScalar hx, PetscScalar hy, PetscScalar hz, int frontfront,
251  |  int frontup, int frontleft, int front, int frontright, int frontdown,
252  |  int upup, int upleft, int up, int upright, int leftleft, int left,
253  |  int current, int right, int rightright, int downleft, int down, int downright,
254  |  int downdown,
255  |  int backup, int backleft, int back, int backright, int backdown,
256  |  int backback)
257  | {
258  |   PetscScalar retval, hx2, hy2, hz2, hxhy, hxhz, hyhz;
259  | 
260  |   hx2 = hx*hx;  hy2 = hy*hy;  hz2 = hz*hz;
261  |   hxhy = hx*hy; hxhz = hx*hz; hyhz = hy*hz;
262  | 
263  |   /*+ This calculates the
264  |     +latex+$\beta$-term, $\kappa\beta$ times the Laplacian of $\Psi'(C)$,
265  |     -latex-beta term, kappa*beta times the laplacian of Psi prime(C),
266  |     +*/
267  |   retval = beta *
268  |     (hx * (ch_psiprime(conc[left], mparam) + ch_psiprime(conc[right], mparam))
269  |      + hy * (ch_psiprime(conc[up], mparam) + ch_psiprime(conc[down], mparam))
270  |      + hz * (ch_psiprime(conc[front], mparam)+ ch_psiprime(conc[back], mparam))
271  |      - 2*(hx+hy+hz) * ch_psiprime(conc[current], mparam));
272  | 
273  |   /*+ then subtracts the
274  |     +latex+$\alpha$-term, $\kappa\alpha\nabla^2\nabla^2C$.
275  |     -latex-alpha term, kappa*alpha times the Laplacian of the Laplacian of C.
276  |     +*/
277  |   retval -= alpha *
278  |     (hx2 * (conc[leftleft] + conc[rightright])
279  |      + hy2 * (conc[upup] + conc[downdown])
280  |      + hz2 * (conc[frontfront] + conc[backback])
281  |      + 2*hxhy * (conc[upleft]+conc[upright]+conc[downleft]+conc[downright])
282  |      + 2*hxhz*(conc[frontleft]+conc[frontright]+conc[backleft]+conc[backright])
283  |      + 2*hyhz * (conc[frontup]+conc[frontdown]+conc[backup]+conc[backdown])
284  |      - (4*hx2 + 4*hxhy + 4*hxhz) * (conc[left] + conc[right])
285  |      - (4*hy2 + 4*hxhy + 4*hyhz) * (conc[up] + conc[down])
286  |      - (4*hz2 + 4*hxhz + 4*hyhz) * (conc[front] + conc[back])
287  |      + (6*hx2 + 6*hy2 + 6*hz2 + 8*hxhy + 8*hxhz + 8*hyhz) * conc[current]);
288  | 
289  |   return retval;
290  | }
291  | 
292  | #define RESIDUAL_FLOPS_3D (7*PSIPRIME_FLOPS + 14 + 65)
293  | 
294  | 
295  | /*++++++++++++++++++++++++++++++++++++++
296  |   This evaluates the nonlinear finite difference approximation to the residuals
297  |   +latex+$R_i$.
298  |   +html+ <i>R<sub>i</sub></i>.
299  |   Note that it loops on the interior points and the boundary separately, to
300  |   avoid conditional statements within the double loop over the local grid
301  |   indices.
302  | 
303  |   int ch_residual_vector_2d It returns zero or an error value.
304  | 
305  |   Vec X The pre-allocated local vector of unknowns.
306  | 
307  |   Vec F The pre-allocated local vector of residuals, filled by this function.
308  | 
309  |   void *ptr Data passed in the application context.
310  |   ++++++++++++++++++++++++++++++++++++++*/
311  | 
312  | int ch_residual_vector_2d (Vec X, Vec F, void *ptr)
313  | {
314  |   AppCtx  *user = (AppCtx*)ptr;
315  |   int     ierr,i,j, mc,chvar, mx,my,xs,ys,xm,ym,gxs,gys,gxm,gym;
316  |   int     xints,xinte,yints,yinte;
317  |   PetscScalar  hx,hy,dhx,dhy,hxm2,hym2;
318  |   PetscScalar  alpha, beta, mparam;
319  |   PetscScalar  *x,*f;
320  |   Vec     localX = user->localX,localF = user->localF; 
321  | 
322  |   /* Scatter ghost points to local vector, using the 2-step process
323  |         DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
324  |      By placing code between these two statements, computations can be
325  |      done while messages are in transition. */
326  |   ierr = DAGlobalToLocalBegin (user->da, X, INSERT_VALUES, localX);
327  |   CHKERRQ (ierr);
328  | 
329  |   /* Define mesh intervals ratios for uniform grid.
330  |      [Note: FD formulae below may someday be normalized by multiplying through
331  |      by local volume element to obtain coefficients O(1) in two dimensions.] */
332  | 
333  |   mc = user->mc;
334  |   chvar = user->chvar;
335  |   mx = user->mx;
336  |   my = user->my;
337  |   mparam = user->mparam;
338  |   alpha = user->kappa * user->gamma * user->epsilon;
339  |   beta = user->kappa * user->gamma / user->epsilon;
340  |   dhx = (PetscScalar)(mx-1);
341  |   dhy = (PetscScalar)(my-1);
342  |   /* These next two lines hard-code the 1x1 square */
343  |   hx = 1./dhx;
344  |   hy = 1./dhy;
345  |   hxm2 = 1./hx/hx;
346  |   hym2 = 1./hy/hy;
347  | 
348  |   ierr = DAGlobalToLocalEnd (user->da, X, INSERT_VALUES, localX);
349  |   CHKERRQ (ierr);
350  | 
351  |   /* Get pointers to vector data */
352  |   ierr = VecGetArray (localX, &x); CHKERRQ (ierr);
353  |   ierr = VecGetArray (localF, &f); CHKERRQ (ierr);
354  | 
355  |   /* Get local grid boundaries */
356  |   ierr = DAGetCorners (user->da, &xs, &ys, PETSC_NULL, &xm, &ym, PETSC_NULL);
357  |   CHKERRQ (ierr);
358  |   ierr = DAGetGhostCorners (user->da, &gxs, &gys, PETSC_NULL, &gxm, &gym,
359  | 			    PETSC_NULL); CHKERRQ (ierr);
360  | 
361  |   /* Determine the interior of the locally owned part of the grid. */
362  |   if (user->period == DA_XPERIODIC || user->period == DA_XYPERIODIC) {
363  |     xints = xs; xinte = xs+xm; }
364  |   else {
365  |     xints = (xs==0) ? 2:xs; xinte = (xs+xm==mx) ? xs+xm-2 : xs+xm; }
366  |   if (user->period == DA_YPERIODIC || user->period == DA_XYPERIODIC) {
367  |     yints = ys; yinte = ys+ym; }
368  |   else {
369  |     yints = (ys==0) ? 2:ys; yinte = (ys+ym==my) ? ys+ym-2 : ys+ym; }
370  | 
371  |   /* Loop over the interior points */
372  |   for (j=yints-gys; j<yinte-gys; j++)
373  |     for (i=j*gxm + xints-gxs; i<j*gxm + xinte-gxs; i++)
374  |       f[C(i)] = ch_residual_2d
375  | 	(x, alpha, beta, mparam, hxm2, hym2,
376  | 	 C(i+2*gxm), C(i+gxm-1), C(i+gxm), C(i+gxm+1),
377  | 	 C(i-2), C(i-1), C(i), C(i+1), C(i+2),
378  | 	 C(i-gxm-1), C(i-gxm), C(i-gxm+1), C(i-2*gxm));
379  | 
380  |   /* Okay, that was the easy part.  Now for the fun part.
381  |      The following conditionals implement symmetry boundary conditions. */
382  | 
383  |   /* Test whether we are on the bottom edge of the global array */
384  |   if (yints == 2) /* Not ys==0 because that would be true for periodic too */
385  |     {
386  |       printf ("This must only be reached if we are NOT y-periodic\n");
387  |       /* Bottom two edge lines */
388  |       for (i=xints-gxs; i<xinte-gxs; i++)
389  | 	{
390  | 	  f[C(i)] = ch_residual_2d
391  | 	    (x, alpha, beta, mparam, hxm2, hym2,
392  | 	     C(i+2*gxm), C(i+gxm-1), C(i+gxm), C(i+gxm+1),
393  | 	     C(i-2), C(i-1), C(i), C(i+1), C(i+2),
394  | 	     C(i+gxm-1), C(i+gxm), C(i+gxm+1), C(i+2*gxm));
395  | 	  f[C(i+gxm)] = ch_residual_2d
396  | 	    (x, alpha, beta, mparam, hxm2, hym2,
397  | 	     C(i+3*gxm), C(i+2*gxm-1), C(i+2*gxm), C(i+2*gxm+1),
398  | 	     C(i+gxm-2), C(i+gxm-1), C(i+gxm), C(i+gxm+1), C(i+gxm+2),
399  | 	     C(i-1), C(i), C(i+1), C(i+gxm));
400  | 	}
401  | 
402  |       /* Bottom left corner */
403  |       if (xints == 2)
404  | 	{
405  | 	  f[C(0)] = ch_residual_2d
406  | 	    (x, alpha, beta, mparam, hxm2, hym2,
407  | 	     C(2*gxm), C(gxm+1), C(gxm), C(gxm+1),
408  | 	     C(2), C(1), C(0), C(1), C(2),
409  | 	     C(gxm+1), C(gxm), C(gxm+1), C(2*gxm));
410  | 	  f[C(1)] = ch_residual_2d
411  | 	    (x, alpha, beta, mparam, hxm2, hym2,
412  | 	     C(2*gxm+1), C(gxm), C(gxm+1), C(gxm+2),
413  | 	     C(1), C(0), C(1), C(2), C(3),
414  | 	     C(gxm), C(gxm+1), C(gxm+2), C(2*gxm+1));
415  | 	  f[C(gxm)] = ch_residual_2d
416  | 	    (x, alpha, beta, mparam, hxm2, hym2,
417  | 	     C(3*gxm), C(2*gxm+1), C(2*gxm), C(2*gxm+1),
418  | 	     C(gxm+2), C(gxm+1), C(gxm+0), C(gxm+1), C(gxm+2),
419  | 	     C(1), C(0), C(1), C(gxm));
420  | 	  f[C(gxm+1)] = ch_residual_2d
421  | 	    (x, alpha, beta, mparam, hxm2, hym2,
422  | 	     C(3*gxm+1), C(2*gxm), C(2*gxm+1), C(2*gxm+2),
423  | 	     C(gxm+1), C(gxm), C(gxm+1), C(gxm+2), C(gxm+3),
424  | 	     C(0), C(1), C(2), C(gxm+1));
425  | 	}
426  | 
427  |       /* Bottom right corner */
428  |       if (xinte == mx-2)
429  | 	{
430  | 	  i = mx-gxs-1; /* Array index of the bottom right corner point */
431  | 	  f[C(i)] = ch_residual_2d
432  | 	    (x, alpha, beta, mparam, hxm2, hym2,
433  | 	     C(i+2*gxm), C(i+gxm-1), C(i+gxm), C(i+gxm-1),
434  | 	     C(i-2), C(i-1), C(i), C(i-1), C(i-2),
435  | 	     C(i+gxm-1), C(i+gxm), C(i+gxm-1), C(i+2*gxm));
436  | 	  f[C(i-1)] = ch_residual_2d
437  | 	    (x, alpha, beta, mparam, hxm2, hym2,
438  | 	     C(i+2*gxm-1), C(i+gxm-2), C(i+gxm-1), C(i+gxm),
439  | 	     C(i-3), C(i-2), C(i-1), C(i), C(i-1),
440  | 	     C(i+gxm-2), C(i+gxm-1), C(i+gxm), C(i+2*gxm-1));
441  | 	  f[C(i+gxm)] = ch_residual_2d
442  | 	    (x, alpha, beta, mparam, hxm2, hym2,
443  | 	     C(i+3*gxm), C(i+2*gxm-1), C(i+2*gxm), C(i+2*gxm-1),
444  | 	     C(i+gxm-2), C(i+gxm-1), C(i+gxm), C(i+gxm-1), C(i+gxm-2),
445  | 	     C(i-1), C(i), C(i-1), C(i+gxm));
446  | 	  f[C(i+gxm-1)] = ch_residual_2d
447  | 	    (x, alpha, beta, mparam, hxm2, hym2,
448  | 	     C(i+3*gxm-1), C(i+2*gxm-2), C(i+2*gxm-1), C(i+2*gxm),
449  | 	     C(i+gxm-3), C(i+gxm-2), C(i+gxm-1), C(i+gxm), C(i+gxm-1),
450  | 	     C(i-2), C(i-1), C(i), C(i+gxm-1));
451  | 	}
452  |     }
453  | 
454  |   /* Test whether we are on the top edge of the global array */
455  |   if (yinte == my-2)
456  |     {
457  |       printf ("This must only be reached if we are NOT y-periodic\n");
458  |       /* Top two edge lines */
459  |       for (i=(my-gys-1)*gxm + xints-gxs; i<(my-gys-1)*gxm + xinte-gxs; i++)
460  | 	{
461  | 	  f[C(i)] = ch_residual_2d
462  | 	    (x, alpha, beta, mparam, hxm2, hym2,
463  | 	     C(i-2*gxm), C(i-gxm-1), C(i-gxm), C(i-gxm+1),
464  | 	     C(i-2), C(i-1), C(i), C(i+1), C(i+2),
465  | 	     C(i-gxm-1), C(i-gxm), C(i-gxm+1), C(i-2*gxm));
466  | 	  f[C(i-gxm)] = ch_residual_2d
467  | 	    (x, alpha, beta, mparam, hxm2, hym2,
468  | 	     C(i-gxm), C(i-1), C(i), C(i+1),
469  | 	     C(i-gxm-2), C(i-gxm-1), C(i-gxm), C(i-gxm+1), C(i-gxm+2),
470  | 	     C(i-2*gxm-1), C(i-2*gxm), C(i-2*gxm+1), C(i-3*gxm));
471  | 	}
472  | 
473  |       /* Top left corner */
474  |       if (xints == 2)
475  | 	{
476  | 	  i = (my-gys-1)*gxm; /* Array index of the top left corner point */
477  | 	  f[C(i)] = ch_residual_2d
478  | 	    (x, alpha, beta, mparam, hxm2, hym2,
479  | 	     C(i-2*gxm), C(i-gxm+1), C(i-gxm), C(i-gxm+1),
480  | 	     C(i+2), C(i+1), C(i), C(i+1), C(i+2),
481  | 	     C(i-gxm+1), C(i-gxm), C(i-gxm+1), C(i-2*gxm));
482  | 	  f[C(i+1)] = ch_residual_2d
483  | 	    (x, alpha, beta, mparam, hxm2, hym2,
484  | 	     C(i-2*gxm+1), C(i-gxm), C(i-gxm+1), C(i-gxm+2),
485  | 	     C(i+1), C(i), C(i+1), C(i+2), C(i+3),
486  | 	     C(i-gxm), C(i-gxm+1), C(i-gxm+2), C(i-2*gxm+1));
487  | 	  f[C(i-gxm)] = ch_residual_2d
488  | 	    (x, alpha, beta, mparam, hxm2, hym2,
489  | 	     C(i-gxm), C(i+1), C(i), C(i+1),
490  | 	     C(i-gxm+2), C(i-gxm+1), C(i-gxm), C(i-gxm+1), C(i-gxm+2),
491  | 	     C(i-2*gxm+1), C(i-2*gxm), C(i-2*gxm+1), C(i-3*gxm));
492  | 	  f[C(i-gxm+1)] = ch_residual_2d
493  | 	    (x, alpha, beta, mparam, hxm2, hym2,
494  | 	     C(i-gxm+1), C(i), C(i+1), C(i+2),
495  | 	     C(i-gxm+1), C(i-gxm), C(i-gxm+1), C(i-gxm+2), C(i-gxm+3),
496  | 	     C(i-2*gxm), C(i-2*gxm+1), C(i-2*gxm+2), C(i-3*gxm+1));
497  | 	}
498  | 
499  |       /* Top right corner */
500  |       if (xinte == mx-2)
501  | 	{ /* Array index of top right corner point */
502  | 	  i = (my-gys-1)*gxm + mx-gxs-1;
503  | 	  f[C(i)] = ch_residual_2d
504  | 	    (x, alpha, beta, mparam, hxm2, hym2,
505  | 	     C(i-2*gxm), C(i-gxm-1), C(i-gxm), C(i-gxm-1),
506  | 	     C(i-2), C(i-1), C(i), C(i-1), C(i-2),
507  | 	     C(i-gxm-1), C(i-gxm), C(i-gxm-1), C(i-2*gxm));
508  | 	  f[C(i-1)] = ch_residual_2d
509  | 	    (x, alpha, beta, mparam, hxm2, hym2,
510  | 	     C(i-2*gxm-1), C(i-gxm-2), C(i-gxm-1), C(i-gxm),
511  | 	     C(i-3), C(i-2), C(i-1), C(i), C(i-1),
512  | 	     C(i-gxm-2), C(i-gxm-1), C(i-gxm), C(i-2*gxm-1));
513  | 	  f[C(i-gxm)] = ch_residual_2d
514  | 	    (x, alpha, beta, mparam, hxm2, hym2,
515  | 	     C(i-gxm), C(i-1), C(i), C(i-1),
516  | 	     C(i-gxm-2), C(i-gxm-1), C(i-gxm), C(i-gxm-1), C(i-gxm-2),
517  | 	     C(i-2*gxm-1), C(i-2*gxm), C(i-2*gxm-1), C(i-3*gxm));
518  | 	  f[C(i-gxm-1)] = ch_residual_2d
519  | 	    (x, alpha, beta, mparam, hxm2, hym2,
520  | 	     C(i-gxm-1), C(i-2), C(i-1), C(i),
521  | 	     C(i-gxm-3), C(i-gxm-2), C(i-gxm-1), C(i-gxm), C(i-gxm-1),
522  | 	     C(i-2*gxm-2), C(i-2*gxm-1), C(i-2*gxm), C(i-3*gxm-1));
523  | 	}
524  |     }
525  | 
526  |   /* Test whether we are on the left edge of the global array */
527  |   if (xints == 2)
528  |     {
529  |       printf ("This must only be reached if we are NOT x-periodic\n");
530  |       /* Left 2 edge lines */
531  |       for (j=yints-gys; j<yinte-gys; j++)
532  | 	{
533  | 	  i = j*gxm;
534  | 	  f[C(i)] = ch_residual_2d
535  | 	    (x, alpha, beta, mparam, hxm2, hym2,
536  | 	     C(i+2*gxm), C(i+gxm+1), C(i+gxm), C(i+gxm+1),
537  | 	     C(i+2), C(i+1), C(i), C(i+1), C(i+2),
538  | 	     C(i-gxm+1), C(i-gxm), C(i-gxm+1), C(i-2*gxm));
539  | 	  f[C(i+1)] = ch_residual_2d
540  | 	    (x, alpha, beta, mparam, hxm2, hym2,
541  | 	     C(i+2*gxm+1), C(i+gxm), C(i+gxm+1), C(i+gxm+2),
542  | 	     C(i+1), C(i), C(i+1), C(i+2), C(i+3),
543  | 	     C(i-gxm), C(i-gxm+1), C(i-gxm+2), C(i-2*gxm+1));
544  | 	}
545  |     }
546  | 
547  |   /* Test whether we are on the right edge of the global array */
548  |   if (xinte == mx-2)
549  |     {
550  |       printf ("This must only be reached if we are NOT x-periodic\n");
551  |       /* Right 2 edge lines */ 
552  |       for (j=yints-gys; j<yinte-gys; j++)
553  | 	{
554  | 	  i = j*gxm + mx-gxs-1;
555  | 	  f[C(i)] = ch_residual_2d
556  | 	    (x, alpha, beta, mparam, hxm2, hym2,
557  | 	     C(i+2*gxm), C(i+gxm-1), C(i+gxm), C(i+gxm-1),
558  | 	     C(i-2), C(i-1), C(i), C(i-1), C(i-2),
559  | 	     C(i-gxm-1), C(i-gxm), C(i-gxm-1), C(i-2*gxm));
560  | 	  f[C(i-1)] = ch_residual_2d
561  | 	    (x, alpha, beta, mparam, hxm2, hym2,
562  | 	     C(i+2*gxm-1), C(i+gxm-2), C(i+gxm-1), C(i+gxm),
563  | 	     C(i-3), C(i-2), C(i-1), C(i), C(i-1),
564  | 	     C(i-gxm-2), C(i-gxm-1), C(i-gxm), C(i-2*gxm-1));
565  | 	}
566  |     }
567  | 
568  |   /* Restore vectors */
569  |   ierr = VecRestoreArray (localX, &x); CHKERRQ (ierr);
570  |   ierr = VecRestoreArray (localF, &f); CHKERRQ (ierr);
571  | 
572  |   /* Insert values into global vector */
573  |   ierr = DALocalToGlobal(user->da,localF,INSERT_VALUES,F); CHKERRQ (ierr);
574  | 
575  |   /* Flop count (multiply-adds are counted as 2 operations) */
576  |   ierr = PetscLogFlops(RESIDUAL_FLOPS_2D*ym*xm); CHKERRQ (ierr);
577  | 
578  |   return 0; 
579  | }
580  | 
581  | 
582  | /*++++++++++++++++++++++++++++++++++++++
583  |   This evaluates the nonlinear finite difference approximation to the residuals
584  |   +latex+$R_i$.
585  |   +html+ <i>R<sub>i</sub></i>.
586  |   Note that it loops on the interior points and the boundary separately, to
587  |   avoid conditional statements within the double loop over the local grid
588  |   indices.
589  |   +latex+\par
590  |   +html+ <p>
591  |   Thus far, only periodic boundary conditions work.
592  | 
593  |   int ch_residual_vector_3d It returns zero or an error value.
594  | 
595  |   Vec X The pre-allocated local vector of unknowns.
596  | 
597  |   Vec F The pre-allocated local vector of residuals, filled by this function.
598  | 
599  |   void *ptr Data passed in the application context.
600  |   ++++++++++++++++++++++++++++++++++++++*/
601  | 
602  | int ch_residual_vector_3d (Vec X, Vec F, void *ptr)
603  | {
604  |   AppCtx  *user = (AppCtx*)ptr;
605  |   int     ierr,i,j,k, mc,chvar, mx,my,mz,xs,ys,zs,xm,ym,zm;
606  |   int     gxs,gys,gzs,gxm,gym,gzm, xints,xinte,yints,yinte,zints,zinte;
607  |   PetscScalar  hx,hy,hz,dhx,dhy,dhz,hxm2,hym2,hzm2;
608  |   PetscScalar  alpha, beta, mparam;
609  |   PetscScalar  *x,*f;
610  |   Vec     localX = user->localX,localF = user->localF; 
611  | 
612  |   /* Scatter ghost points to local vector, using the 2-step process
613  |         DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
614  |      By placing code between these two statements, computations can be
615  |      done while messages are in transition. */
616  |   ierr = DAGlobalToLocalBegin (user->da, X, INSERT_VALUES, localX);
617  |   CHKERRQ (ierr);
618  | 
619  |   /* Define mesh intervals ratios for uniform grid.
620  |      [Note: FD formulae below may someday be normalized by multiplying through
621  |      by local volume element to obtain coefficients O(1) in two dimensions.] */
622  | 
623  |   mc = user->mc;
624  |   chvar = user->chvar;
625  |   mx = user->mx; my = user->my; mz = user->mz;
626  |   mparam = user->mparam;
627  |   alpha = user->kappa * user->gamma * user->epsilon;
628  |   beta = user->kappa * user->gamma / user->epsilon;
629  |   dhx = (PetscScalar)(mx-1); dhy = (PetscScalar)(my-1);
630  |   dhz = (PetscScalar)(mz-1);
631  |   /* This next line hard-codes the 1x1x1 cube */
632  |   hx = 1./dhx;     hy = 1./dhy;     hz = 1./dhz;
633  |   hxm2 = 1./hx/hx; hym2 = 1./hy/hy; hzm2 = 1./hz/hz;
634  | 
635  |   ierr = DAGlobalToLocalEnd (user->da, X, INSERT_VALUES, localX);
636  |   CHKERRQ (ierr);
637  | 
638  |   /* Get pointers to vector data */
639  |   ierr = VecGetArray (localX, &x); CHKERRQ (ierr);
640  |   ierr = VecGetArray (localF, &f); CHKERRQ (ierr);
641  | 
642  |   /* Get local grid boundaries */
643  |   ierr = DAGetCorners (user->da, &xs, &ys, &zs, &xm, &ym, &zm); CHKERRQ (ierr);
644  |   ierr = DAGetGhostCorners (user->da, &gxs, &gys, &gzs, &gxm, &gym, &gzm);
645  |   CHKERRQ (ierr);
646  | 
647  |   /* Determine the interior of the locally owned part of the grid. */
648  |   if (user->period == DA_XYZPERIODIC) {
649  |     xints = xs; xinte = xs+xm;
650  |     yints = ys; yinte = ys+ym;
651  |     zints = zs; zinte = zs+zm; }
652  |   else {
653  |     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,
654  | 	    "Only periodic boundary conditions in 3-D Cahn-Hilliard"); }
655  | 
656  |   /* Loop over the interior points (no boundaries yet) */
657  |   for (k=zints-gzs; k<zinte-gzs; k++)
658  |     for (j=k*gym + yints-gys; j<k*gym + yinte-gys; j++)
659  |       for (i=j*gxm + xints-gxs; i<j*gxm + xinte-gxs; i++)
660  | 	f[C(i)] = ch_residual_3d
661  | 	  (x, alpha, beta, mparam, hxm2, hym2, hzm2,
662  | 	   C(i+2*gxm*gym), C(i+gxm*gym+gxm), C(i+gxm*gym-1), C(i+gxm*gym),
663  | 	   C(i+gxm*gym+1), C(i+gxm*gym-gxm),
664  | 	   C(i+2*gxm), C(i+gxm-1), C(i+gxm), C(i+gxm+1),
665  | 	   C(i-2), C(i-1), C(i), C(i+1), C(i+2),
666  | 	   C(i-gxm-1), C(i-gxm), C(i-gxm+1), C(i-2*gxm),
667  | 	   C(i-gxm*gym+gxm), C(i-gxm*gym-1), C(i-gxm*gym), C(i-gxm*gym+1),
668  | 	   C(i-gxm*gym-gxm), C(i-2*gxm*gym));
669  | 
670  |   /* Restore vectors */
671  |   ierr = VecRestoreArray (localX, &x); CHKERRQ (ierr);
672  |   ierr = VecRestoreArray (localF, &f); CHKERRQ (ierr);
673  | 
674  |   /* Insert values into global vector */
675  |   ierr = DALocalToGlobal(user->da,localF,INSERT_VALUES,F); CHKERRQ (ierr);
676  | 
677  |   /* Flop count (multiply-adds are counted as 2 operations) */
678  |   ierr = PetscLogFlops(RESIDUAL_FLOPS_3D*ym*xm*zm); CHKERRQ (ierr);
679  | 
680  |   /* ierr = VecView (F, VIEWER_STDOUT_SELF); CHKERRQ (ierr); */
681  | 
682  |   return 0;
683  | }
684  | 
685  | 
686  | /*++++++++++++++++++++++++++++++++++++++
687  |   This computes the Jacobian matrix at each iteration, starting with the alpha
688  |   term which is pre-computed at the beginning by
689  |   +latex+{\tt ch\_jacobian\_alpha\_2d()}.
690  |   +html+ <tt>ch_jacobian_alpha_2d()</tt>.
691  | 
692  |   int ch_jacobian_2d It returns 0 or an error code.
693  | 
694  |   Vec X The vector of unknowns.
695  | 
696  |   Mat *A The Jacobian matrix returned to PETSc.
697  | 
698  |   Mat *B The matrix preconditioner, in this case the same matrix.
699  | 
700  |   MatStructure *flag Flag saying the nonzeroes are the same each time.
701  | 
702  |   void *ptr Application context structure.
703  |   ++++++++++++++++++++++++++++++++++++++*/
704  | 
705  | int ch_jacobian_2d (Vec X, Mat *A, Mat *B, MatStructure *flag, void *ptr)
706  | {
707  |   AppCtx  *user = (AppCtx*)ptr;
708  |   int     ierr,node,i,j, mc,chvar, mx,my,xs,ys,xm,ym,gxs,gys,gxm,gym;
709  |   int     xints,xinte,yints,yinte;
710  |   int     columns [5];
711  |   PetscScalar  hx,hy,dhx,dhy,hxm2,hym2;
712  |   PetscScalar  alpha,beta,mparam;
713  |   PetscScalar  *x, bvalue[5];
714  |   Vec     localX = user->localX;
715  | 
716  |   /* Set the MatStructure flag right, Mats to return */
717  |   *flag = SAME_NONZERO_PATTERN;
718  |   A = &(user->J);
719  |   B = &(user->J);
720  | 
721  |   /* Copy the alpha term Jacobian into the main Jacobian space */
722  |   ierr = MatCopy (user->alpha, user->J, SAME_NONZERO_PATTERN);
723  | 
724  |   /* Scatter ghost points to local vector, using 2-step async I/O with (a
725  |      couple of trivial) computations in between. */
726  | 
727  |   ierr = DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX); CHKERRQ (ierr);
728  |   mc = user->mc;  chvar = user->chvar;
729  |   mx = user->mx;  my = user->my;  mparam = user->mparam;
730  |   alpha = user->kappa * user->gamma * user->epsilon;
731  |   beta = user->kappa * user->gamma / user->epsilon;
732  | 
733  |   /* Define mesh intervals ratios for uniform grid.
734  |      [Note: FD formulae below may someday be normalized by multiplying through
735  |      by local volume element to obtain coefficients O(1) in two dimensions.] */
736  |   dhx = (PetscScalar)(mx-1);  dhy = (PetscScalar)(my-1);
737  |   hx = 1./dhx;           hy = 1./dhy; /* This line hard-codes the 1x1 square */
738  |   hxm2 = 1./hx/hx;       hym2 = 1./hy/hy;
739  | 
740  |   ierr = DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX); CHKERRQ (ierr);
741  | 
742  |   /* Get pointer to vector data */
743  |   ierr = VecGetArray(localX,&x); CHKERRQ (ierr);
744  | 
745  |   /*
746  |      Get local grid boundaries
747  |   */
748  |   ierr = DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ (ierr);
749  |   ierr = DAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); CHKERRQ (ierr);
750  | 
751  |   /* Determine the interior of the locally owned part of the grid. */
752  |   if (user->period == DA_XPERIODIC || user->period == DA_XYPERIODIC) {
753  |     xints = xs; xinte = xs+xm; }
754  |   else {
755  |     xints = (xs==0) ? 1:xs; xinte = (xs+xm==mx) ? xs+xm-1 : xs+xm; }
756  |   if (user->period == DA_YPERIODIC || user->period == DA_XYPERIODIC) {
757  |     yints = ys; yinte = ys+ym; }
758  |   else {
759  |     yints = (ys==0) ? 1:ys; yinte = (ys+ym==my) ? ys+ym-1 : ys+ym; }
760  | 
761  |   /* Loop over the interior points... */
762  |   for (j=yints-gys; j<yinte-gys; j++)
763  |     for (i=j*gxm + xints-gxs; i<j*gxm + xinte-gxs; i++) {
764  | 
765  |       /* Form the column indices */
766  |       columns[0]=C(i+gxm);
767  |       columns[1]=C(i-1); columns[2]=C(i); columns[3]=C(i+1);
768  |       columns[4]=C(i-gxm);
769  | 
770  |       bvalue[0] = beta * hym2 * ch_psidoubleprime (x[columns[0]], mparam);
771  |       bvalue[1] = beta * hxm2 * ch_psidoubleprime (x[columns[1]], mparam);
772  |       bvalue[2] = -2.*beta * (hxm2+hym2) *
773  | 	ch_psidoubleprime (x[columns[2]], mparam);
774  |       bvalue[3] = beta * hxm2 * ch_psidoubleprime (x[columns[3]], mparam);
775  |       bvalue[4] = beta * hym2 * ch_psidoubleprime (x[columns[4]], mparam);
776  | 
777  |       node = C(i);
778  |       MatSetValuesLocal (user->J, 1,&node, 5,columns, bvalue, ADD_VALUES);
779  |     }
780  | 
781  |   /* Now the boundary conditions... */
782  | 
783  |   /* Test whether we're on the bottom row and non-y-periodic */
784  |   if (yints == 1) {
785  |     printf ("This must only be reached if we are NOT y-periodic\n");
786  |     for (i=xints-gxs; i<xinte-gxs; i++) {
787  | 
788  |       /* Form the column indices */
789  |       columns[0]=C(i+gxm);
790  |       columns[1]=C(i-1); columns[2]=C(i); columns[3]=C(i+1);
791  | 
792  |       bvalue[0] = 2.*beta * hym2 * ch_psidoubleprime (x[columns[0]], mparam);
793  |       bvalue[1] = beta * hxm2 * ch_psidoubleprime (x[columns[1]], mparam);
794  |       bvalue[2] = -2.*beta * (hxm2+hym2) *
795  | 	ch_psidoubleprime (x[columns[2]], mparam);
796  |       bvalue[3] = beta * hxm2 * ch_psidoubleprime (x[columns[3]], mparam);
797  | 
798  |       node = C(i);
799  |       MatSetValuesLocal (user->J, 1,&node, 4,columns, bvalue, ADD_VALUES);
800  |     }
801  | 
802  |     /* Bottom left corner */
803  |     if (xints == 1) {
804  |       columns[0]=C(i=0); columns[1]=C(i+1);
805  |       columns[2]=C(i+gxm);
806  | 
807  |       bvalue[0] = -2.*beta * (hxm2+hym2) *
808  | 	ch_psidoubleprime (x[columns[0]], mparam);
809  |       bvalue[1] = 2.*beta * hxm2 * ch_psidoubleprime (x[columns[1]], mparam);
810  |       bvalue[2] = 2.*beta * hym2 * ch_psidoubleprime (x[columns[2]], mparam);
811  | 
812  |       node = C(i);
813  |       MatSetValuesLocal (user->J, 1,&node, 3,columns, bvalue, ADD_VALUES);
814  |     }
815  | 
816  |     /* Bottom right corner */
817  |     if (xinte == mx-1) {
818  |       columns[0]=C(i=mx-gxs-1); columns[1]=C(i-1);
819  |       columns[2]=C(i+gxm);
820  | 
821  |       bvalue[0] = -2.*beta * (hxm2+hym2) *
822  | 	ch_psidoubleprime (x[columns[0]], mparam);
823  |       bvalue[1] = 2.*beta * hxm2 * ch_psidoubleprime (x[columns[1]], mparam);
824  |       bvalue[2] = 2.*beta * hym2 * ch_psidoubleprime (x[columns[2]], mparam);
825  | 
826  |       node = C(i);
827  |       MatSetValuesLocal (user->J, 1,&node, 3,columns, bvalue, ADD_VALUES);
828  |     }
829  |   }
830  | 
831  |   /* Test whether we're on the top row and non-y-periodic */
832  |   if (yinte == my-1) {
833  |     printf ("This must only be reached if we are NOT y-periodic\n");
834  |     for (i=(my-gys-1)*gxm + xints-gxs; i<(my-gys-1)*gxm + xinte-gxs; i++) {
835  | 
836  |       /* Form the column indices */
837  |       columns[0]=C(i-1); columns[1]=C(i); columns[2]=C(i+1);
838  |       columns[3]=C(i-gxm);
839  | 
840  |       bvalue[0] = beta * hxm2 * ch_psidoubleprime (x[columns[0]], mparam);
841  |       bvalue[1] = -2.*beta * (hxm2+hym2) *
842  | 	ch_psidoubleprime (x[columns[1]], mparam);
843  |       bvalue[2] = beta * hxm2 * ch_psidoubleprime (x[columns[2]], mparam);
844  |       bvalue[3] = 2.*beta * hym2 * ch_psidoubleprime (x[columns[3]], mparam);
845  | 
846  |       node = C(i);
847  |       MatSetValuesLocal (user->J, 1,&node, 4,columns, bvalue, ADD_VALUES);
848  |     }
849  | 
850  |     /* Top left corner */
851  |     if (xints == 1) {
852  |       columns[0]=C(i=(my-gys-1)*gxm); columns[1] = C(i+1);
853  |       columns[2]=C(i-gxm);
854  | 
855  |       bvalue[0] = -2.*beta * (hxm2+hym2) *
856  | 	ch_psidoubleprime (x[columns[0]], mparam);
857  |       bvalue[1] = 2.*beta * hxm2 * ch_psidoubleprime (x[columns[1]], mparam);
858  |       bvalue[2] = 2.*beta * hym2 * ch_psidoubleprime (x[columns[2]], mparam);
859  | 
860  |       node = C(i);
861  |       MatSetValuesLocal (user->J, 1,&node, 3,columns, bvalue, ADD_VALUES);
862  |     }
863  | 
864  |     /* Top right corner */
865  |     if (xinte == mx-1) {
866  |       columns[0]=C(i=(my-gys-1)*gxm + mx-gxs-1); columns[1]=C(i-1);
867  |       columns[2]=C(i-gxm);
868  | 
869  |       bvalue[0] = -2.*beta * (hxm2+hym2) *
870  | 	ch_psidoubleprime (x[columns[0]], mparam);
871  |       bvalue[1] = 2.*beta * hxm2 * ch_psidoubleprime (x[columns[1]], mparam);
872  |       bvalue[2] = 2.*beta * hym2 * ch_psidoubleprime (x[columns[2]], mparam);
873  | 
874  |       node = C(i);
875  |       MatSetValuesLocal (user->J, 1,&node, 3,columns, bvalue, ADD_VALUES);
876  |     }
877  |   }
878  | 
879  |   /* Left column */
880  |   if (xints == 1) {
881  |     printf ("This must only be reached if we are NOT x-periodic\n");
882  |     for (i=(yints-gys)*gxm; i<(yinte-gys)*gxm; i+=gxm) {
883  |       columns[0]=C(i+gxm);
884  |       columns[1]=C(i); columns[2] = C(i+1);
885  |       columns[3]=C(i-gxm);
886  | 
887  |       bvalue[0] = beta * hym2 * ch_psidoubleprime (x[columns[0]], mparam);
888  |       bvalue[1] = -2.*beta * (hxm2+hym2) *
889  | 	ch_psidoubleprime (x[columns[1]], mparam);
890  |       bvalue[2] = 2.*beta * hxm2 * ch_psidoubleprime (x[columns[2]], mparam);
891  |       bvalue[3] = beta * hym2 * ch_psidoubleprime (x[columns[3]], mparam);
892  | 
893  |       node = C(i);
894  |       MatSetValuesLocal (user->J, 1,&node, 4,columns, bvalue, ADD_VALUES);
895  |     }
896  |   }
897  | 
898  |   /* Right column */
899  |   if (xinte == mx-1) {
900  |     printf ("This must only be reached if we are NOT x-periodic\n");
901  |     for (i=(yints-gys)*gxm + mx-gxs-1; i<(yinte-gys)*gxm + mx-gxs-1; i+=gxm) {
902  |       columns[0]=C(i+gxm);
903  |       columns[1]=C(i-1); columns[2] = C(i);
904  |       columns[3]=C(i-gxm);
905  | 
906  |       bvalue[0] = beta * hym2 * ch_psidoubleprime (x[columns[0]], mparam);
907  |       bvalue[1] = 2.*beta * hxm2 * ch_psidoubleprime (x[columns[1]], mparam);
908  |       bvalue[2] = -2.*beta * (hxm2+hym2) *
909  | 	ch_psidoubleprime (x[columns[2]], mparam);
910  |       bvalue[3] = beta * hym2 * ch_psidoubleprime (x[columns[3]], mparam);
911  | 
912  |       node = C(i);
913  |       MatSetValuesLocal (user->J, 1,&node, 4,columns, bvalue, ADD_VALUES);
914  |     }
915  |   }
916  | 
917  |   /* Assemble matrix, using the 2-step process:
918  |      MatAssemblyBegin(), MatAssemblyEnd(). */
919  |   ierr = MatAssemblyBegin(user->J,MAT_FINAL_ASSEMBLY); CHKERRQ (ierr);
920  |   ierr = MatAssemblyEnd(user->J,MAT_FINAL_ASSEMBLY); CHKERRQ (ierr);
921  | 
922  |   /* Restore unknown vector */
923  |   ierr = VecRestoreArray(localX,&x); CHKERRQ (ierr);
924  | 
925  |   return ierr;
926  | }
927  | 
928  | 
929  | /*++++++++++++++++++++++++++++++++++++++
930  |   This computes the Jacobian matrix at each iteration, starting with the alpha
931  |   term which is pre-computed at the beginning by
932  |   +latex+{\tt ch\_jacobian\_alpha\_3d()}.
933  |   +html+ <tt>ch_jacobian_alpha_3d()</tt>.
934  | 
935  |   int ch_jacobian_3d It returns 0 or an error code.
936  | 
937  |   Vec X The vector of unknowns.
938  | 
939  |   Mat *A The Jacobian matrix returned to PETSc.
940  | 
941  |   Mat *B The matrix preconditioner, in this case the same matrix.
942  | 
943  |   MatStructure *flag Flag saying the nonzeroes are the same each time.
944  | 
945  |   void *ptr Application context structure.
946  |   ++++++++++++++++++++++++++++++++++++++*/
947  | 
948  | int ch_jacobian_3d (Vec X, Mat *A, Mat *B, MatStructure *flag, void *ptr)
949  | {
950  |   AppCtx  *user = (AppCtx*)ptr;
951  |   int     ierr,node,i,j,k, mc,chvar, mx,my,mz, xs,ys,zs, xm,ym,zm;
952  |   int     gxs,gys,gzs, gxm,gym,gzm, xints,xinte, yints,yinte, zints,zinte;
953  |   int     columns [7];
954  |   PetscScalar  hx,hy,hz, dhx,dhy,dhz, hxm2,hym2,hzm2;
955  |   PetscScalar  alpha,beta,mparam;
956  |   PetscScalar  *x, bvalue[7];
957  |   Vec     localX = user->localX;
958  | 
959  |   /* Set the MatStructure flag right, Mats to return */
960  |   *flag = SAME_NONZERO_PATTERN;
961  |   A = &(user->J);
962  |   B = &(user->J);
963  | 
964  |   /* Copy the alpha term Jacobian into the main Jacobian space */
965  |   ierr = MatCopy (user->alpha, user->J, SAME_NONZERO_PATTERN);
966  | 
967  |   /* Scatter ghost points to local vector, using 2-step async I/O with (a
968  |      couple of trivial) computations in between. */
969  | 
970  |   ierr = DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX); CHKERRQ (ierr);
971  |   mc = user->mc;  chvar = user->chvar;
972  |   mx = user->mx;  my = user->my;  mz = user->mz;  mparam = user->mparam;
973  |   alpha = user->kappa * user->gamma * user->epsilon;
974  |   beta = user->kappa * user->gamma / user->epsilon;
975  | 
976  |   /* Define mesh intervals ratios for uniform grid.
977  |      [Note: FD formulae below may someday be normalized by multiplying through
978  |      by local volume element to obtain coefficients O(1) in two dimensions.] */
979  |   dhx = (PetscScalar)(mx-1);  dhy = (PetscScalar)(my-1);
980  |   dhz = (PetscScalar)(mz-1);
981  |   /* This line hard-codes the 1x1x1 cube */
982  |   hx = 1./dhx;           hy = 1./dhy;          hz = 1./dhz;
983  |   hxm2 = 1./hx/hx;       hym2 = 1./hy/hy;      hzm2 = 1./hz/hz;
984  | 
985  |   ierr = DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX); CHKERRQ (ierr);
986  | 
987  |   /* Get pointer to vector data */
988  |   ierr = VecGetArray(localX,&x); CHKERRQ (ierr);
989  | 
990  |   /* Get local grid boundaries */
991  |   ierr = DAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm); CHKERRQ (ierr);
992  |   ierr = DAGetGhostCorners(user->da,&gxs,&gys,&gzs,&gxm,&gym,&gzm);
993  |   CHKERRQ (ierr);
994  | 
995  |   /* Determine the interior of the locally owned part of the grid. */
996  |   if (user->period == DA_XYZPERIODIC) {
997  |     xints = xs; xinte = xs+xm;
998  |     yints = ys; yinte = ys+ym;
999  |     zints = zs; zinte = zs+zm; }
1000 |   else {
1001 |     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,
1002 | 	    "Only periodic boundary conditions in 3-D Cahn-Hilliard"); }
1003 | 
1004 |   /* Loop over the interior points... */
1005 |   for (k=zints-gzs; k<zinte-gzs; k++)
1006 |     for (j=k*gym + yints-gys; j<k*gym + yinte-gys; j++)
1007 |       for (i=j*gxm + xints-gxs; i<j*gxm + xinte-gxs; i++) {
1008 | 
1009 | 	/* Form the column indices */
1010 | 	columns[0]=C(i+gxm*gym); columns[1]=C(i+gxm);
1011 | 	columns[2]=C(i-1); columns[3]=C(i); columns[4]=C(i+1);
1012 | 	columns[5]=C(i-gxm); columns[6]=C(i-gxm*gym);
1013 | 
1014 | 	bvalue[0] = beta * hzm2 * ch_psidoubleprime (x[columns[0]], mparam);
1015 | 	bvalue[1] = beta * hym2 * ch_psidoubleprime (x[columns[1]], mparam);
1016 | 	bvalue[2] = beta * hxm2 * ch_psidoubleprime (x[columns[2]], mparam);
1017 | 	bvalue[3] = -2.*beta * (hxm2+hym2+hzm2) *
1018 | 	  ch_psidoubleprime (x[columns[3]], mparam);
1019 | 	bvalue[4] = beta * hxm2 * ch_psidoubleprime (x[columns[4]], mparam);
1020 | 	bvalue[5] = beta * hym2 * ch_psidoubleprime (x[columns[5]], mparam);
1021 | 	bvalue[6] = beta * hzm2 * ch_psidoubleprime (x[columns[6]], mparam);
1022 | 
1023 | 	node = C(i);
1024 | 	MatSetValuesLocal (user->J, 1,&node, 5,columns, bvalue, ADD_VALUES);
1025 |       }
1026 | 
1027 |   /* Assemble matrix, using the 2-step process:
1028 |      MatAssemblyBegin(), MatAssemblyEnd(). */
1029 |   ierr = MatAssemblyBegin (user->J,MAT_FINAL_ASSEMBLY); CHKERRQ (ierr);
1030 |   ierr = MatAssemblyEnd (user->J,MAT_FINAL_ASSEMBLY); CHKERRQ (ierr);
1031 | 
1032 |   /* Restore unknown vector */
1033 |   ierr = VecRestoreArray (localX,&x); CHKERRQ (ierr);
1034 | 
1035 |   /* ierr = MatView (user->J, VIEWER_STDOUT_SELF); CHKERRQ (ierr); */
1036 | 
1037 |   return ierr;
1038 | }
1039 | 
1040 | 
1041 | /*++++++++++++++++++++++++++++++++++++++
1042 |   This creates the initial alpha and J matrices, where alpha is the alpha term
1043 |   component of the Jacobian.  Since the alpha term is linear, this part of the
1044 |   Jacobian need only be calculated once.
1045 | 
1046 |   int ch_jacobian_alpha_2d It returns zero or an error message.
1047 | 
1048 |   AppCtx *user The application context structure pointer.
1049 |   ++++++++++++++++++++++++++++++++++++++*/
1050 | 
1051 | int ch_jacobian_alpha_2d (AppCtx *user)
1052 | {
1053 |   int     ierr,node,i,j, mc,chvar, mx,my,xs,ys,xm,ym,gxs,gys,gxm,gym;
1054 |   int     xints,xinte,yints,yinte;
1055 |   int     columns [13];
1056 |   PetscScalar  hx,hy,dhx,dhy,hxm2,hym2;
1057 |   PetscScalar  alpha,beta,mparam;
1058 |   PetscScalar  avalue [25];
1059 | 
1060 |   mc = user->mc;  chvar = user->chvar;
1061 |   mx = user->mx;  my = user->my;  mparam = user->mparam;
1062 |   alpha = user->kappa * user->gamma * user->epsilon;
1063 |   beta = user->kappa * user->gamma / user->epsilon;
1064 | 
1065 |   /* Define mesh intervals ratios for uniform grid.
1066 |      [Note: FD formulae below may someday be normalized by multiplying through
1067 |      by local volume element to obtain coefficients O(1) in two dimensions.] */
1068 |   dhx = (PetscScalar)(mx-1);  dhy = (PetscScalar)(my-1);
1069 |   hx = 1./dhx;           hy = 1./dhy; /* This line hard-codes the 1x1 square */
1070 |   hxm2 = 1./hx/hx;       hym2 = 1./hy/hy;
1071 | 
1072 |   /* Get local grid boundaries */
1073 |   ierr = DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
1074 |   CHKERRQ (ierr);
1075 |   ierr = DAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL);
1076 |   CHKERRQ (ierr);
1077 | 
1078 |   /* Determine the interior of the locally owned part of the grid. */
1079 |   if (user->period == DA_XPERIODIC || user->period == DA_XYPERIODIC) {
1080 |     xints = xs; xinte = xs+xm; }
1081 |   else {
1082 |     xints = (xs==0) ? 2:xs; xinte = (xs+xm==mx) ? xs+xm-2 : xs+xm; }
1083 |   if (user->period == DA_YPERIODIC || user->period == DA_XYPERIODIC) {
1084 |     yints = ys; yinte = ys+ym; }
1085 |   else {
1086 |     yints = (ys==0) ? 2:ys; yinte = (ys+ym==my) ? ys+ym-2 : ys+ym; }
1087 | 
1088 |   /* Get parameters for matrix creation */
1089 |   i = xm * ym * user->mc;
1090 |   j = mx * my * user->mc;
1091 | 
1092 |   /* Create the distributed alpha matrix */
1093 |   ierr = MatCreateMPIAIJ (user->comm, i,i, j,j, 13,PETSC_NULL, 6,PETSC_NULL,
1094 | 			  &(user->alpha));
1095 |   CHKERRQ (ierr);
1096 | 
1097 |   /* Get the local-to-global mapping and associate it with the matrix */
1098 |   ierr = DAGetISLocalToGlobalMapping (user->da, &user->isltog); CHKERRQ (ierr);
1099 |   ierr = MatSetLocalToGlobalMapping (user->alpha, user->isltog); CHKERRQ (ierr);
1100 | 
1101 |   /* Pre-compute alpha term values (see ch_residual_2d above)
1102 |      This should be the only place they're actually computed; they'll be
1103 |      used below. */
1104 |   avalue[0] = avalue[12] = -alpha*hym2*hym2;
1105 |   avalue[1] = avalue[3] = avalue[9] = avalue[11] = -2.*alpha*hxm2*hym2;
1106 |   avalue[2] = avalue[10] = 4.*alpha*hym2*(hxm2+hym2);
1107 |   avalue[4] = avalue[8] = -alpha*hxm2*hxm2;
1108 |   avalue[5] = avalue[7] = 4.*alpha*hxm2*(hxm2+hym2);
1109 |   avalue[6] = -alpha*(6.*hxm2*hxm2 + 6.*hym2*hym2 + 8.*hxm2*hym2);
1110 | 
1111 |   /* Loop over interior nodes */
1112 |   for (j=yints-gys; j<yinte-gys; j++)
1113 |     for (i=j*gxm + xints-gxs; i<j*gxm + xinte-gxs; i++) {
1114 | 
1115 |       /* Form the column indices */
1116 |       columns[0]=C(i+2*gxm);
1117 |       columns[1]=C(i+gxm-1); columns[2]=C(i+gxm);  columns[3]=C(i+gxm+1);
1118 |       columns[4]=C(i-2);     columns[5]=C(i-1);    columns[6]=C(i);
1119 |       columns[7]=C(i+1);     columns[8]=C(i+2);
1120 |       columns[9]=C(i-gxm-1); columns[10]=C(i-gxm); columns[11]=C(i-gxm+1);
1121 |       columns[12]=C(i-2*gxm);
1122 | 
1123 |       node = C(i);
1124 |       MatSetValuesLocal (user->alpha, 1,&node, 13,columns, avalue,
1125 | 			 INSERT_VALUES);
1126 |     }
1127 | 
1128 |   /* Okay, that was the easy part, now for the boundary conditions! */
1129 | 
1130 |   /* Test whether we're on the bottom row and non-y-periodic */
1131 |   if (yints == 2) {
1132 |     printf ("This must only be reached if we are NOT y-periodic\n");
1133 |     /* Change value 6 and do the second-from-bottom row */
1134 |     avalue[6] += avalue[12];
1135 |     for (i=gxm + xints-gxs; i<gxm + xinte-gxs; i++) {
1136 | 
1137 |       /* Form the column indices */
1138 |       columns[0]=C(i+2*gxm);
1139 |       columns[1]=C(i+gxm-1); columns[2]=C(i+gxm); columns[3]=C(i+gxm+1);
1140 |       columns[4]=C(i-2);     columns[5]=C(i-1);   columns[6]=C(i);
1141 |       columns[7]=C(i+1);     columns[8]=C(i+2);
1142 |       columns[9]=C(i-gxm-1); columns[10]=C(i-gxm); columns[11]=C(i-gxm+1);
1143 | 
1144 |       node = C(i);
1145 |       MatSetValuesLocal (user->alpha, 1,&node, 12,columns, avalue,
1146 | 			 INSERT_VALUES);
1147 |     }
1148 |     avalue[6] -= avalue[12];
1149 | 
1150 |     /* Make some new avalues and do the bottom row */
1151 |     avalue[13] = 2.*avalue[0];
1152 |     avalue[14] = avalue[16] = 2.*avalue[1];
1153 |     avalue[15] = 2.*avalue[2];
1154 |     avalue[17] = avalue[21] = avalue[4];
1155 |     avalue[18] = avalue[20] = avalue[5];
1156 |     avalue[19] = avalue[6];
1157 |     for (i=xints-gxs; i<xinte-gxs; i++) {
1158 | 
1159 |       /* Form the column indices */
1160 |       columns[0]=C(i+2*gxm);
1161 |       columns[1]=C(i+gxm-1); columns[2]=C(i+gxm); columns[3]=C(i+gxm+1);
1162 |       columns[4]=C(i-2);     columns[5]=C(i-1);   columns[6]=C(i);
1163 |       columns[7]=C(i+1);     columns[8]=C(i+2);
1164 | 
1165 |       node = C(i);
1166 |       MatSetValuesLocal (user->alpha, 1,&node, 9,columns, avalue+13,
1167 | 			 INSERT_VALUES);
1168 |     }
1169 | 
1170 |     /* Bottom left corner */
1171 |     if (xints == 2) {
1172 |       node=C(0); /* The point (0,0) */
1173 |       columns[0]=C(2*gxm);
1174 |       columns[1]=C(gxm);   columns[2]=C(gxm+1);
1175 |       columns[3]=C(0);     columns[4]=C(1);     columns[5]=C(2);
1176 |       avalue[13] = 2.*avalue[0];
1177 |       avalue[14] = 2.*avalue[2];
1178 |       avalue[15] = 4.*avalue[3];
1179 |       avalue[16] = avalue[6];
1180 |       avalue[17] = 2.*avalue[7];
1181 |       avalue[18] = 2.*avalue[8];
1182 |       MatSetValuesLocal (user->alpha, 1,&node, 6,columns, avalue+13,
1183 | 			 INSERT_VALUES);
1184 |       node=C(1); /* The point (1,0) */
1185 |       columns[0]=C(2*gxm+1);
1186 |       columns[1]=C(gxm); columns[2]=C(gxm+1); columns[3]=C(gxm+2);
1187 |       columns[4]=C(0);   columns[5]=C(1);  columns[6]=C(2);  columns[7]=C(3);
1188 |       avalue[13] = 2.*avalue[0];
1189 |       avalue[14] = avalue[16] = 2.*avalue[1];
1190 |       avalue[15] = 2.*avalue[2];
1191 |       avalue[17] = avalue[19] = avalue[5];
1192 |       avalue[18] = avalue[6]  + avalue[4];
1193 |       avalue[20] = avalue[8];
1194 |       MatSetValuesLocal (user->alpha, 1,&node, 8,columns, avalue+13,
1195 | 			 INSERT_VALUES);
1196 |       node=C(gxm); /* The point (0,1) */
1197 |       columns[0]=C(3*gxm);
1198 |       columns[1]=C(2*gxm); columns[2]=C(2*gxm+1);
1199 |       columns[3]=C(gxm);   columns[4]=C(gxm+1);   columns[5]=C(gxm+2);
1200 |       columns[6]=C(0);     columns[7]=C(1);
1201 |       avalue[13] = avalue[0];
1202 |       avalue[14] = avalue[19] = avalue[2];
1203 |       avalue[15] = avalue[20] = 2.*avalue[3];
1204 |       avalue[16] = avalue[6] + avalue[12];
1205 |       avalue[17] = 2.*avalue[7];
1206 |       avalue[18] = 2.*avalue[8];
1207 |       MatSetValuesLocal (user->alpha, 1,&node, 8,columns, avalue+13,
1208 | 			 INSERT_VALUES);
1209 |       node=C(gxm+1); /* The point (1,1) */
1210 |       columns[0]=C(3*gxm+1);
1211 |       columns[1]=C(2*gxm); columns[2]=C(2*gxm+1); columns[3]=C(2*gxm+2);
1212 |       columns[4]=C(gxm);   columns[5]=C(gxm+1);   columns[6]=C(gxm+2);
1213 |       columns[7]=C(gxm+3);
1214 |       columns[8]=C(0);     columns[9]=C(1);       columns[10]=C(2);
1215 |       avalue[13] = avalue[0];
1216 |       avalue[14] = avalue[16] = avalue[21] = avalue[23] = avalue[1];
1217 |       avalue[15] = avalue[22] = avalue[2];
1218 |       avalue[17] = avalue[19] = avalue[5];
1219 |       avalue[18] = avalue[6] + avalue[4] + avalue[12];
1220 |       avalue[20] = avalue[8];
1221 |       MatSetValuesLocal (user->alpha, 1,&node, 11,columns, avalue+13,
1222 | 			 INSERT_VALUES);
1223 |     }
1224 | 
1225 |     /* Bottom right corner */
1226 |     if(xinte == mx-2) {
1227 |       node=C(i=mx-gxs-1); /* The point (mx-1, 0) */
1228 |       columns[0]=C(i+2*gxm);
1229 |       columns[1]=C(i+gxm-1); columns[2]=C(i+gxm);
1230 |       columns[3]=C(i-2);     columns[4]=C(i-1);   columns[5]=C(i);
1231 |       avalue[13] = 2.*avalue[0];
1232 |       avalue[14] = 4.*avalue[1];
1233 |       avalue[15] = 2.*avalue[2];
1234 |       avalue[16] = 2.*avalue[4];
1235 |       avalue[17] = 2.*avalue[5];
1236 |       avalue[18] = avalue[6];
1237 |       MatSetValuesLocal (user->alpha, 1,&node, 6,columns, avalue+13,
1238 | 			 INSERT_VALUES);
1239 |       node=C(i=mx-gxs-2); /* The point (mx-2,0) */
1240 |       columns[0]=C(i+2*gxm);
1241 |       columns[1]=C(i+gxm-1); columns[2]=C(i+gxm); columns[3]=C(i+gxm+1);
1242 |       columns[4]=C(i-2);     columns[5]=C(i-1);   columns[6]=C(i);
1243 |       columns[7]=C(i+1);
1244 |       avalue[13] = 2.*avalue[0];
1245 |       avalue[14] = avalue[16] = 2.*avalue[1];
1246 |       avalue[15] = 2.*avalue[2];
1247 |       avalue[17] = avalue[4];
1248 |       avalue[18] = avalue[20] = avalue[5];
1249 |       avalue[19] = avalue[6]  + avalue[8];
1250 |       MatSetValuesLocal (user->alpha, 1,&node, 8,columns, avalue+13,
1251 | 			 INSERT_VALUES);
1252 |       node=C(i=gxm+mx-gxs-1); /* The point (mx-1,1) */
1253 |       columns[0]=C(i+2*gxm);
1254 |       columns[1]=C(i+gxm-1); columns[2]=C(i+gxm);
1255 |       columns[3]=C(i-2);     columns[4]=C(i-1);   columns[5]=C(i);
1256 |       columns[6]=C(i-gxm-1); columns[7]=C(i-gxm);
1257 |       avalue[13] = avalue[0];
1258 |       avalue[14] = avalue[19] = 2.*avalue[1];
1259 |       avalue[15] = avalue[20] = avalue[2];
1260 |       avalue[16] = 2.*avalue[4];
1261 |       avalue[17] = 2.*avalue[5];
1262 |       avalue[18] = avalue[6] + avalue[12];
1263 |       MatSetValuesLocal (user->alpha, 1,&node, 8,columns, avalue+13,
1264 | 			 INSERT_VALUES);
1265 |       node=C(i=gxm+mx-gxs-2); /* The point (mx-2,1) */
1266 |       columns[0]=C(i+2*gxm);
1267 |       columns[1]=C(i+gxm-1); columns[2]=C(i+gxm); columns[3]=C(i+gxm+1);
1268 |       columns[4]=C(i-2);     columns[5]=C(i-1);   columns[6]=C(i);
1269 |       columns[7]=C(i+1);
1270 |       columns[8]=C(i-gxm-1); columns[9]=C(i-gxm);   columns[10]=C(i-gxm+1);
1271 |       avalue[13] = avalue[0];
1272 |       avalue[14] = avalue[16] = avalue[21] = avalue[23] = avalue[1];
1273 |       avalue[15] = avalue[22] = avalue[2];
1274 |       avalue[17] = avalue[4];
1275 |       avalue[18] = avalue[20] = avalue[5];
1276 |       avalue[19] = avalue[6] + avalue[8] + avalue[12];
1277 |       MatSetValuesLocal (user->alpha, 1,&node, 11,columns, avalue+13,
1278 | 			 INSERT_VALUES);
1279 |     }
1280 |   }
1281 | 
1282 |   /* Test whether we're on the top row and non-y-periodic */
1283 |   if (yinte == my-2) {
1284 |     printf ("This must only be reached if we are NOT y-periodic\n");
1285 |     /* Change value 6 and do the second-from-top row */
1286 |     avalue[6] += avalue[0];
1287 |     for (i=(my-gys-2)*gxm + xints-gxs;
1288 | 	 i<(my-gys-2)*gxm + xinte-gxs; i++) {
1289 | 
1290 |       /* Form the column indices */
1291 |       columns[1]=C(i+gxm-1); columns[2]=C(i+gxm);  columns[3]=C(i+gxm+1);
1292 |       columns[4]=C(i-2);     columns[5]=C(i-1);    columns[6]=C(i);
1293 |       columns[7]=C(i+1);     columns[8]=C(i+2);
1294 |       columns[9]=C(i-gxm-1); columns[10]=C(i-gxm); columns[11]=C(i-gxm+1);
1295 |       columns[12]=C(i-2*gxm);
1296 | 
1297 |       node = C(i);
1298 |       MatSetValuesLocal (user->alpha, 1,&node, 12,columns+1, avalue+1,
1299 | 			 INSERT_VALUES);
1300 |     }
1301 |     avalue[6] -= avalue[0];
1302 | 
1303 |     /* Make some new avalues and do the top row */
1304 |     avalue[13] = avalue[17] = avalue[4];
1305 |     avalue[14] = avalue[16] = avalue[5];
1306 |     avalue[15] = avalue[6];
1307 |     avalue[18] = avalue[20] = 2.*avalue[9];
1308 |     avalue[19] = 2.*avalue[10];
1309 |     avalue[21] = 2.*avalue[12];
1310 |     for (i=(my-gys-1)*gxm + xints-gxs;
1311 | 	 i<(my-gys-1)*gxm + xinte-gxs; i++) {
1312 | 
1313 |       /* Form the column indices */
1314 |       columns[0]=C(i-2);     columns[1]=C(i-1);   columns[2]=C(i);
1315 |       columns[3]=C(i+1);     columns[4]=C(i+2);
1316 |       columns[5]=C(i-gxm-1); columns[6]=C(i-gxm); columns[7]=C(i-gxm+1);
1317 |       columns[8]=C(i-2*gxm);
1318 | 
1319 |       node = C(i);
1320 |       MatSetValuesLocal (user->alpha, 1,&node, 9,columns, avalue+13,
1321 | 			 INSERT_VALUES);
1322 |     }
1323 | 
1324 |     /* Top left corner */
1325 |     if (xints == 2) {
1326 |       node=C(i=(my-gys-1)*gxm); /* The point (0,my-1) */
1327 |       columns[0]=C(i);       columns[1]=C(i+1); columns[2]=C(i+2);
1328 |       columns[3]=C(i-gxm);   columns[4]=C(i-gxm+1);
1329 |       columns[5]=C(i-2*gxm);
1330 |       avalue[13] = avalue[6];
1331 |       avalue[14] = 2.*avalue[7];
1332 |       avalue[15] = 2.*avalue[8];
1333 |       avalue[16] = 2.*avalue[10];
1334 |       avalue[17] = 4.*avalue[11];
1335 |       avalue[18] = 2.*avalue[12];
1336 |       MatSetValuesLocal (user->alpha, 1,&node, 6,columns, avalue+13,
1337 | 			 INSERT_VALUES);
1338 |       node=C(i=(my-gys-1)*gxm+1); /* The point (1,my-1) */
1339 |       columns[0]=C(i-1);     columns[1]=C(i);     columns[2]=C(i+1);
1340 |       columns[3]=C(i+2);
1341 |       columns[4]=C(i-gxm-1); columns[5]=C(i-gxm); columns[6]=C(i-gxm+1);
1342 |       columns[7]=C(i-2*gxm);
1343 |       avalue[13] = avalue[15] = avalue[5];
1344 |       avalue[14] = avalue[6]  + avalue[4];
1345 |       avalue[16] = avalue[8];
1346 |       avalue[17] = avalue[19] = 2.*avalue[9];
1347 |       avalue[18] = 2.*avalue[10];
1348 |       avalue[20] = 2.*avalue[12];
1349 |       MatSetValuesLocal (user->alpha, 1,&node, 8,columns, avalue+13,
1350 | 			 INSERT_VALUES);
1351 |       node=C(i=(my-gys-2)*gxm); /* The point (0,my-2) */
1352 |       columns[0]=C(i+gxm); columns[1]=C(i+gxm+1);
1353 |       columns[2]=C(i);     columns[3]=C(i+1);   columns[4]=C(i+2);
1354 |       columns[5]=C(i-gxm); columns[6]=C(i-gxm+1);
1355 |       columns[7]=C(i-2*gxm);
1356 |       avalue[13] = avalue[18] = avalue[2];
1357 |       avalue[14] = avalue[19] = 2.*avalue[3];
1358 |       avalue[15] = avalue[6] + avalue[0];
1359 |       avalue[16] = 2.*avalue[7];
1360 |       avalue[17] = 2.*avalue[8];
1361 |       avalue[20] = avalue[12];
1362 |       MatSetValuesLocal (user->alpha, 1,&node, 8,columns, avalue+13,
1363 | 			 INSERT_VALUES);
1364 |       node=C(i=(my-gys-2)*gxm+1); /* The point (1,my-2) */
1365 |       columns[0]=C(i+gxm-1); columns[1]=C(i+gxm); columns[2]=C(i+gxm+1);
1366 |       columns[3]=C(i-1);     columns[4]=C(i);     columns[5]=C(i+1);
1367 |       columns[6]=C(i+2);
1368 |       columns[7]=C(i-gxm-1); columns[8]=C(i-gxm); columns[9]=C(i-gxm+1);
1369 |       columns[10]=C(i-2*gxm);
1370 |       avalue[13] = avalue[15] = avalue[20] = avalue[22] = avalue[1];
1371 |       avalue[14] = avalue[21] = avalue[2];
1372 |       avalue[16] = avalue[18] = avalue[5];
1373 |       avalue[17] = avalue[6] + avalue[4] + avalue[0];
1374 |       avalue[19] = avalue[8];
1375 |       avalue[23] = avalue[12];
1376 |       MatSetValuesLocal (user->alpha, 1,&node, 11,columns, avalue+13,
1377 | 			 INSERT_VALUES);
1378 |     }
1379 | 
1380 |     /* Top right corner */
1381 |     if(xinte == mx-2) {
1382 |       node=C(i=(my-gys-1)*gxm + mx-gxs-1); /* The point (mx-1, my-1) */
1383 |       columns[0]=C(i-2);   columns[1]=C(i-1);   columns[2]=C(i);
1384 |       columns[3]=C(i-gxm-1); columns[4]=C(i-gxm);
1385 |       columns[5]=C(i-2*gxm);
1386 |       avalue[13] = 2.*avalue[4];
1387 |       avalue[14] = 2.*avalue[5];
1388 |       avalue[15] = avalue[6];
1389 |       avalue[16] = 4.*avalue[9];
1390 |       avalue[17] = 2.*avalue[10];
1391 |       avalue[18] = 2.*avalue[12];
1392 |       MatSetValuesLocal (user->alpha, 1,&node, 6,columns, avalue+13,
1393 | 			 INSERT_VALUES);
1394 |       node=C(i=(my-gys-1)*gxm + mx-gxs-2); /* The point (mx-2, my-1) */
1395 |       columns[0]=C(i-2);   columns[1]=C(i-1);   columns[2]=C(i);
1396 |       columns[3]=C(i+1);
1397 |       columns[4]=C(i-gxm-1); columns[5]=C(i-gxm); columns[6]=C(i-gxm+1);
1398 |       columns[7]=C(i-2*gxm);
1399 |       avalue[13] = avalue[4];
1400 |       avalue[14] = avalue[16] = avalue[5];
1401 |       avalue[15] = avalue[6]  + avalue[8];
1402 |       avalue[17] = avalue[19] = 2.*avalue[9];
1403 |       avalue[18] = 2.*avalue[10];
1404 |       avalue[20] = 2.*avalue[12];
1405 |       MatSetValuesLocal (user->alpha, 1,&node, 8,columns, avalue+13,
1406 | 			 INSERT_VALUES);
1407 |       node=C(i=(my-gys-2)*gxm + mx-gxs-1); /* The point (mx-1, my-2) */
1408 |       columns[0]=C(i+gxm-1); columns[1]=C(i+gxm);
1409 |       columns[2]=C(i-2);     columns[3]=C(i-1);   columns[4]=C(i);
1410 |       columns[5]=C(i-gxm-1); columns[6]=C(i-gxm);
1411 |       columns[7]=C(i-2*gxm);
1412 |       avalue[13] = avalue[18] = 2.*avalue[1];
1413 |       avalue[14] = avalue[19] = avalue[2];
1414 |       avalue[15] = 2.*avalue[4];
1415 |       avalue[16] = 2.*avalue[5];
1416 |       avalue[17] = avalue[6] + avalue[0];
1417 |       avalue[20] = avalue[12];
1418 |       MatSetValuesLocal (user->alpha, 1,&node, 8,columns, avalue+13,
1419 | 			 INSERT_VALUES);
1420 |       node=C(i=(my-gys-2)*gxm + mx-gxs-2); /* The point (mx-2, my-2) */
1421 |       columns[1]=C(i+gxm-1); columns[2]=C(i+gxm); columns[3]=C(i+gxm+1);
1422 |       columns[4]=C(i-2);     columns[5]=C(i-1);   columns[6]=C(i);
1423 |       columns[7]=C(i+1);
1424 |       columns[8]=C(i-gxm-1); columns[9]=C(i-gxm); columns[10]=C(i-gxm+1);
1425 |       columns[11]=C(i-2*gxm);
1426 |       avalue[14] = avalue[16] = avalue[21] = avalue[23] = avalue[1];
1427 |       avalue[15] = avalue[22] = avalue[2];
1428 |       avalue[17] = avalue[4];
1429 |       avalue[18] = avalue[20] = avalue[5];
1430 |       avalue[19] = avalue[6] + avalue[8] + avalue[0];
1431 |       avalue[24] = avalue[12];
1432 |       MatSetValuesLocal (user->alpha, 1,&node, 11,columns+1, avalue+14,
1433 | 			 INSERT_VALUES);
1434 |     }
1435 |   }
1436 | 
1437 |   /* Test whether we're on the left column and non-x-periodic */
1438 |   if (xints == 2) {
1439 |     printf ("This must only be reached if we are NOT x-periodic\n");
1440 |     /* Make some avalues and do the second-from-left column */
1441 |     avalue[13] = avalue[24] = avalue[0];
1442 |     avalue[14] = avalue[16] = avalue[21] = avalue[23] = avalue[1];
1443 |     avalue[15] = avalue[22] = avalue[2];
1444 |     avalue[17] = avalue[19] = avalue[5];
1445 |     avalue[18] = avalue[6]  + avalue[4];
1446 |     avalue[20] = avalue[8];
1447 |     for (i=(yints-gys)*gxm + 1; i<(yinte-gys)*gxm + 1; i+=gxm) {
1448 |       columns[0]=C(i+2*gxm);
1449 |       columns[1]=C(i+gxm-1); columns[2]=C(i+gxm);  columns[3]=C(i+gxm+1);
1450 |       columns[4]=C(i-1);     columns[5]=C(i);      columns[6]=C(i+1);
1451 |       columns[7]=C(i+2);
1452 |       columns[8]=C(i-gxm-1); columns[9]=C(i-gxm); columns[10]=C(i-gxm+1);
1453 |       columns[11]=C(i-2*gxm);
1454 | 
1455 |       node = C(i);
1456 |       MatSetValuesLocal (user->alpha, 1,&node, 12,columns, avalue+13,
1457 | 			 INSERT_VALUES);
1458 |     }
1459 | 
1460 |     /* Make some more avalues and do the leftmost column */
1461 |     avalue[13] = avalue[21] = avalue[0];
1462 |     avalue[14] = avalue[19] = avalue[2];
1463 |     avalue[15] = avalue[20] = 2.*avalue[3];
1464 |     avalue[16] = avalue[6];
1465 |     avalue[17] = 2.*avalue[7];
1466 |     avalue[18] = 2.*avalue[8];
1467 |     for (i=(yints-gys)*gxm; i<(yinte-gys)*gxm; i+=gxm) {
1468 |       columns[0]=C(i+2*gxm);
1469 |       columns[1]=C(i+gxm);  columns[2]=C(i+gxm+1);
1470 |       columns[3]=C(i);      columns[4]=C(i+1);     columns[5]=C(i+2);
1471 |       columns[6]=C(i-gxm);  columns[7]=C(i-gxm+1);
1472 |       columns[8]=C(i-2*gxm);
1473 | 
1474 |       node = C(i);
1475 |       MatSetValuesLocal (user->alpha, 1,&node, 9,columns, avalue+13,
1476 | 			 INSERT_VALUES);
1477 |     }
1478 |   }
1479 | 
1480 |   /* Test whether we're on the right column and non-x-periodic */
1481 |   if (xinte == mx-2) {
1482 |     printf ("This must only be reached if we are NOT x-periodic\n");
1483 |     /* Make some avalues and do the second-from-right column */
1484 |     avalue[13] = avalue[24] = avalue[0];
1485 |     avalue[14] = avalue[16] = avalue[21] = avalue[23] = avalue[1];
1486 |     avalue[15] = avalue[22] = avalue[2];
1487 |     avalue[17] = avalue[4];
1488 |     avalue[18] = avalue[20] = avalue[5];
1489 |     avalue[19] = avalue[6]  + avalue[8];
1490 |     for (i=(yints-gys)*gxm + mx-gxs-2; i<(yinte-gys)*gxm + mx-gxs-2; i+=gxm) {
1491 |       columns[0]=C(i+2*gxm);
1492 |       columns[1]=C(i+gxm-1); columns[2]=C(i+gxm);  columns[3]=C(i+gxm+1);
1493 |       columns[4]=C(i-2);     columns[5]=C(i-1);    columns[6]=C(i);
1494 |       columns[7]=C(i+1);
1495 |       columns[8]=C(i-gxm-1); columns[9]=C(i-gxm); columns[10]=C(i-gxm+1);
1496 |       columns[11]=C(i-2*gxm);
1497 | 
1498 |       node = C(i);
1499 |       MatSetValuesLocal (user->alpha, 1,&node, 12,columns, avalue+13,
1500 | 			 INSERT_VALUES);
1501 |     }
1502 | 
1503 |     /* Make some more avalues and do the rightmost column */
1504 |     avalue[13] = avalue[21] = avalue[0];
1505 |     avalue[14] = avalue[19] = 2.*avalue[1];
1506 |     avalue[15] = avalue[20] = avalue[2];
1507 |     avalue[16] = 2.*avalue[4];
1508 |     avalue[17] = 2.*avalue[5];
1509 |     avalue[18] = avalue[6];
1510 |     for (i=(yints-gys)*gxm + mx-gxs-1; i<(yinte-gys)*gxm + mx-gxs-1; i+=gxm) {
1511 |       columns[0]=C(i+2*gxm);
1512 |       columns[1]=C(i+gxm-1);  columns[2]=C(i+gxm);
1513 |       columns[3]=C(i-2);      columns[4]=C(i-1);     columns[5]=C(i);
1514 |       columns[6]=C(i-gxm-1);  columns[7]=C(i-gxm);
1515 |       columns[8]=C(i-2*gxm);
1516 | 
1517 |       node = C(i);
1518 |       MatSetValuesLocal (user->alpha, 1,&node, 9,columns, avalue+13,
1519 | 			 INSERT_VALUES);
1520 |     }
1521 |   }
1522 | 
1523 |   /* Assemble matrix, using the 2-step process:
1524 |      MatAssemblyBegin(), MatAssemblyEnd(). */
1525 |   ierr = MatAssemblyBegin(user->alpha,MAT_FINAL_ASSEMBLY); CHKERRQ (ierr);
1526 |   ierr = MatAssemblyEnd(user->alpha,MAT_FINAL_ASSEMBLY); CHKERRQ (ierr);
1527 | 
1528 |   /* Make J a copy of alpha with the same local mapping (setting new mapping is
1529 |      unnecessary with PETSc 2.1.5). */
1530 |   ierr = MatDuplicate (user->alpha, MAT_COPY_VALUES, &(user->J)); CHKERRQ (ierr);
1531 | 
1532 |   /* Tell the matrix we will never add a new nonzero location to the
1533 |      matrix. If we do, it will generate an error. */
1534 |   ierr = MatSetOption(user->J,MAT_NEW_NONZERO_LOCATION_ERR); CHKERRQ (ierr);
1535 | 
1536 |   return ierr;
1537 | }
1538 | 
1539 | 
1540 | /*++++++++++++++++++++++++++++++++++++++
1541 |   This creates the initial alpha and J matrices, where alpha is the alpha term
1542 |   component of the Jacobian.  Since the alpha term is linear, this part of the
1543 |   Jacobian need only be calculated once.
1544 | 
1545 |   int ch_jacobian_alpha_3d It returns zero or an error message.
1546 | 
1547 |   AppCtx *user The application context structure pointer.
1548 |   ++++++++++++++++++++++++++++++++++++++*/
1549 | 
1550 | int ch_jacobian_alpha_3d (AppCtx *user)
1551 | {
1552 |   int     ierr,node,i,j,k, mc,chvar, mx,my,mz, xs,ys,zs, xm,ym,zm;
1553 |   int     gxs,gys,gzs, gxm,gym,gzm, xints,xinte, yints,yinte, zints,zinte;
1554 |   int     columns [25];
1555 |   PetscScalar  hx,hy,hz, dhx,dhy,dhz, hxm2,hym2,hzm2;
1556 |   PetscScalar  alpha,beta,mparam;
1557 |   PetscScalar  avalue [25];
1558 | 
1559 |   mc = user->mc;  chvar = user->chvar;
1560 |   mx = user->mx;  my = user->my;  mz = user->mz;  mparam = user->mparam;
1561 |   alpha = user->kappa * user->gamma * user->epsilon;
1562 |   beta = user->kappa * user->gamma / user->epsilon;
1563 | 
1564 |   /* Define mesh intervals ratios for uniform grid.
1565 |      [Note: FD formulae below may someday be normalized by multiplying through
1566 |      by local volume element to obtain coefficients O(1) in two dimensions.] */
1567 |   dhx = (PetscScalar)(mx-1);  dhy = (PetscScalar)(my-1);
1568 |   dhz = (PetscScalar)(mz-1);
1569 |   /* This line hard-codes the 1x1x1 cube */
1570 |   hx = 1./dhx;           hy = 1./dhy;          hz = 1./dhz;
1571 |   hxm2 = 1./hx/hx;       hym2 = 1./hy/hy;      hzm2 = 1./hz/hz;
1572 | 
1573 |   /* Get local grid boundaries */
1574 |   ierr = DAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm); CHKERRQ (ierr);
1575 |   ierr = DAGetGhostCorners(user->da,&gxs,&gys,&gzs,&gxm,&gym,&gzm);
1576 |   CHKERRQ (ierr);
1577 | 
1578 |   /* Determine the interior of the locally owned part of the grid. */
1579 |   if (user->period == DA_XYZPERIODIC) {
1580 |     xints = xs; xinte = xs+xm;
1581 |     yints = ys; yinte = ys+ym;
1582 |     zints = zs; zinte = zs+zm; }
1583 |   else {
1584 |     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,
1585 | 	    "Only periodic boundary conditions in 3-D Cahn-Hilliard"); }
1586 | 
1587 |   /* Get parameters for matrix creation */
1588 |   i = xm * ym * zm * user->mc;
1589 |   j = mx * my * mz * user->mc;
1590 | 
1591 |   /* Create the distributed alpha matrix */
1592 |   ierr = MatCreateMPIAIJ (user->comm, i,i, j,j, 25,PETSC_NULL, 12,PETSC_NULL,
1593 | 			  &(user->alpha));
1594 |   CHKERRQ (ierr);
1595 | 
1596 |   /* Get the local-to-global mapping and associate it with the matrix */
1597 |   ierr = DAGetISLocalToGlobalMapping (user->da, &user->isltog); CHKERRQ (ierr);
1598 |   ierr = MatSetLocalToGlobalMapping (user->alpha, user->isltog); CHKERRQ (ierr);
1599 | 
1600 |   /* Pre-compute alpha term values (see ch_residual_2d above)
1601 |      This should be the only place they're actually computed; they'll be
1602 |      used below. */
1603 |   avalue[0]  = avalue[24] = -alpha*hzm2*hzm2;
1604 |   avalue[1]  = avalue[5] = avalue[19] = avalue[23] = -2.*alpha*hym2*hzm2;
1605 |   avalue[2]  = avalue[4] = avalue[20] = avalue[22] = -2.*alpha*hxm2*hzm2;
1606 |   avalue[3]  = avalue[21] = 4.*alpha*hzm2*(hxm2+hym2+hzm2);
1607 |   avalue[6]  = avalue[18] = -alpha*hym2*hym2;
1608 |   avalue[7]  = avalue[9] = avalue[15] = avalue[17] = -2.*alpha*hxm2*hym2;
1609 |   avalue[8]  = avalue[16] = 4.*alpha*hym2*(hxm2+hym2+hzm2);
1610 |   avalue[10] = avalue[14] = -alpha*hxm2*hxm2;
1611 |   avalue[11] = avalue[13] = 4.*alpha*hxm2*(hxm2+hym2+hzm2);
1612 |   avalue[12] = -alpha*(6.*hxm2*hxm2 + 6.*hym2*hym2 + 6.*hzm2*hzm2
1613 | 		       + 8.*hxm2*hym2 + 8.*hxm2*hzm2 + 8.*hym2*hzm2);
1614 | 
1615 |   /* Loop over interior nodes */
1616 |   for (k=zints-gzs; k<zinte-gzs; k++)
1617 |     for (j=k*gym + yints-gys; j<k*gym + yinte-gys; j++)
1618 |       for (i=j*gxm + xints-gxs; i<j*gxm + xinte-gxs; i++) {
1619 | 
1620 | 	/* Form the column indices */
1621 | 	columns[0]=C(i+2*gxm*gym);
1622 | 	columns[1]=C(i+gxm*gym+gxm);
1623 | 	columns[2]=C(i+gxm*gym-1);
1624 | 	columns[3]=C(i+gxm*gym);
1625 | 	columns[4]=C(i+gxm*gym+1);
1626 | 	columns[5]=C(i+gxm*gym-gxm);
1627 | 	columns[6]=C(i+2*gxm);
1628 | 	columns[7]=C(i+gxm-1);  columns[8]=C(i+gxm);  columns[9]=C(i+gxm+1);
1629 | 	columns[10]=C(i-2);     columns[11]=C(i-1);   columns[12]=C(i);
1630 | 	columns[13]=C(i+1);     columns[14]=C(i+2);
1631 | 	columns[15]=C(i-gxm-1); columns[16]=C(i-gxm); columns[17]=C(i-gxm+1);
1632 | 	columns[18]=C(i-2*gxm);
1633 | 	columns[19]=C(i-gxm*gym+gxm);
1634 | 	columns[20]=C(i-gxm*gym-1);
1635 | 	columns[21]=C(i-gxm*gym);
1636 | 	columns[22]=C(i-gxm*gym+1);
1637 | 	columns[23]=C(i-gxm*gym-gxm);
1638 | 	columns[24]=C(i-2*gxm*gym);
1639 | 
1640 | 	node = C(i);
1641 | 	MatSetValuesLocal (user->alpha, 1,&node, 25,columns, avalue,
1642 | 			   INSERT_VALUES);
1643 |       }
1644 | 
1645 |   /* Assemble matrix, using the 2-step process:
1646 |      MatAssemblyBegin(), MatAssemblyEnd(). */
1647 |   ierr = MatAssemblyBegin(user->alpha,MAT_FINAL_ASSEMBLY); CHKERRQ (ierr);
1648 |   ierr = MatAssemblyEnd(user->alpha,MAT_FINAL_ASSEMBLY); CHKERRQ (ierr);
1649 | 
1650 |   /* Make J a copy of alpha with the same local mapping (setting new mapping is
1651 |      unnecessary with PETSc 2.1.5). */
1652 |   ierr = MatDuplicate (user->alpha, MAT_COPY_VALUES, &(user->J)); CHKERRQ (ierr);
1653 | 
1654 |   /* Tell the matrix we will never add a new nonzero location to the
1655 |      matrix. If we do, it will generate an error. */
1656 |   ierr = MatSetOption(user->J,MAT_NEW_NONZERO_LOCATION_ERR); CHKERRQ (ierr);
1657 | 
1658 |   return ierr;
1659 | }