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>>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,PETSC_TRUE); 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,PETSC_TRUE); CHKERRQ (ierr);
1657 |
1658 | return ierr;
1659 | }