Actual source code: ex20.c
2: static char help[] ="Nonlinear Radiative Transport PDE with multigrid in 3d.\n\
3: Uses 3-dimensional distributed arrays.\n\
4: A 3-dim simplified Radiative Transport test problem is used, with analytic Jacobian. \n\
5: \n\
6: Solves the linear systems via multilevel methods \n\
7: \n\
8: The command line\n\
9: options are:\n\
10: -tleft <tl>, where <tl> indicates the left Diriclet BC \n\
11: -tright <tr>, where <tr> indicates the right Diriclet BC \n\
12: -beta <beta>, where <beta> indicates the exponent in T \n\n";
14: /*
16: This example models the partial differential equation
18: - Div(alpha* T^beta (GRAD T)) = 0.
20: where beta = 2.5 and alpha = 1.0
22: BC: T_left = 1.0, T_right = 0.1, dT/dn_top = dTdn_bottom = dT/dn_up = dT/dn_down = 0.
24: in the unit square, which is uniformly discretized in each of x and
25: y in this simple encoding. The degrees of freedom are cell centered.
27: A finite volume approximation with the usual 7-point stencil
28: is used to discretize the boundary value problem to obtain a
29: nonlinear system of equations.
31: This code was contributed by Nickolas Jovanovic based on ex18.c
33: */
35: #include <petscsnes.h>
36: #include <petscdm.h>
37: #include <petscdmda.h>
39: /* User-defined application context */
41: typedef struct {
42: PetscReal tleft,tright; /* Dirichlet boundary conditions */
43: PetscReal beta,bm1,coef; /* nonlinear diffusivity parameterizations */
44: } AppCtx;
46: #define POWFLOP 5 /* assume a pow() takes five flops */
48: extern PetscErrorCode FormInitialGuess(SNES,Vec,void*);
49: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
50: extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
52: int main(int argc,char **argv)
53: {
54: SNES snes;
55: AppCtx user;
56: PetscInt its,lits;
57: PetscReal litspit;
58: DM da;
60: PetscInitialize(&argc,&argv,NULL,help);
61: /* set problem parameters */
62: user.tleft = 1.0;
63: user.tright = 0.1;
64: user.beta = 2.5;
65: PetscOptionsGetReal(NULL,NULL,"-tleft",&user.tleft,NULL);
66: PetscOptionsGetReal(NULL,NULL,"-tright",&user.tright,NULL);
67: PetscOptionsGetReal(NULL,NULL,"-beta",&user.beta,NULL);
68: user.bm1 = user.beta - 1.0;
69: user.coef = user.beta/2.0;
71: /*
72: Set the DMDA (grid structure) for the grids.
73: */
74: DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,5,5,5,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,0,&da);
75: DMSetFromOptions(da);
76: DMSetUp(da);
77: DMSetApplicationContext(da,&user);
79: /*
80: Create the nonlinear solver
81: */
82: SNESCreate(PETSC_COMM_WORLD,&snes);
83: SNESSetDM(snes,da);
84: SNESSetFunction(snes,NULL,FormFunction,&user);
85: SNESSetJacobian(snes,NULL,NULL,FormJacobian,&user);
86: SNESSetFromOptions(snes);
87: SNESSetComputeInitialGuess(snes,FormInitialGuess,NULL);
89: SNESSolve(snes,NULL,NULL);
90: SNESGetIterationNumber(snes,&its);
91: SNESGetLinearSolveIterations(snes,&lits);
92: litspit = ((PetscReal)lits)/((PetscReal)its);
93: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
94: PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %D\n",lits);
95: PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / SNES = %e\n",(double)litspit);
97: SNESDestroy(&snes);
98: DMDestroy(&da);
99: PetscFinalize();
100: return 0;
101: }
102: /* -------------------- Form initial approximation ----------------- */
103: PetscErrorCode FormInitialGuess(SNES snes,Vec X,void *ctx)
104: {
105: AppCtx *user;
106: PetscInt i,j,k,xs,ys,xm,ym,zs,zm;
107: PetscScalar ***x;
108: DM da;
111: SNESGetDM(snes,&da);
112: DMGetApplicationContext(da,&user);
113: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
114: DMDAVecGetArray(da,X,&x);
116: /* Compute initial guess */
117: for (k=zs; k<zs+zm; k++) {
118: for (j=ys; j<ys+ym; j++) {
119: for (i=xs; i<xs+xm; i++) {
120: x[k][j][i] = user->tleft;
121: }
122: }
123: }
124: DMDAVecRestoreArray(da,X,&x);
125: return 0;
126: }
127: /* -------------------- Evaluate Function F(x) --------------------- */
128: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
129: {
130: AppCtx *user = (AppCtx*)ptr;
131: PetscInt i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm;
132: PetscScalar zero = 0.0,one = 1.0;
133: PetscScalar hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy;
134: PetscScalar t0,tn,ts,te,tw,an,as,ae,aw,dn,ds,de,dw,fn = 0.0,fs = 0.0,fe =0.0,fw = 0.0;
135: PetscScalar tleft,tright,beta,td,ad,dd,fd=0.0,tu,au,du=0.0,fu=0.0;
136: PetscScalar ***x,***f;
137: Vec localX;
138: DM da;
141: SNESGetDM(snes,&da);
142: DMGetLocalVector(da,&localX);
143: DMDAGetInfo(da,NULL,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
144: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1); hz = one/(PetscReal)(mz-1);
145: hxhydhz = hx*hy/hz; hyhzdhx = hy*hz/hx; hzhxdhy = hz*hx/hy;
146: tleft = user->tleft; tright = user->tright;
147: beta = user->beta;
149: /* Get ghost points */
150: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
151: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
152: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
153: DMDAVecGetArray(da,localX,&x);
154: DMDAVecGetArray(da,F,&f);
156: /* Evaluate function */
157: for (k=zs; k<zs+zm; k++) {
158: for (j=ys; j<ys+ym; j++) {
159: for (i=xs; i<xs+xm; i++) {
160: t0 = x[k][j][i];
162: if (i > 0 && i < mx-1 && j > 0 && j < my-1 && k > 0 && k < mz-1) {
164: /* general interior volume */
166: tw = x[k][j][i-1];
167: aw = 0.5*(t0 + tw);
168: dw = PetscPowScalar(aw,beta);
169: fw = dw*(t0 - tw);
171: te = x[k][j][i+1];
172: ae = 0.5*(t0 + te);
173: de = PetscPowScalar(ae,beta);
174: fe = de*(te - t0);
176: ts = x[k][j-1][i];
177: as = 0.5*(t0 + ts);
178: ds = PetscPowScalar(as,beta);
179: fs = ds*(t0 - ts);
181: tn = x[k][j+1][i];
182: an = 0.5*(t0 + tn);
183: dn = PetscPowScalar(an,beta);
184: fn = dn*(tn - t0);
186: td = x[k-1][j][i];
187: ad = 0.5*(t0 + td);
188: dd = PetscPowScalar(ad,beta);
189: fd = dd*(t0 - td);
191: tu = x[k+1][j][i];
192: au = 0.5*(t0 + tu);
193: du = PetscPowScalar(au,beta);
194: fu = du*(tu - t0);
196: } else if (i == 0) {
198: /* left-hand (west) boundary */
199: tw = tleft;
200: aw = 0.5*(t0 + tw);
201: dw = PetscPowScalar(aw,beta);
202: fw = dw*(t0 - tw);
204: te = x[k][j][i+1];
205: ae = 0.5*(t0 + te);
206: de = PetscPowScalar(ae,beta);
207: fe = de*(te - t0);
209: if (j > 0) {
210: ts = x[k][j-1][i];
211: as = 0.5*(t0 + ts);
212: ds = PetscPowScalar(as,beta);
213: fs = ds*(t0 - ts);
214: } else {
215: fs = zero;
216: }
218: if (j < my-1) {
219: tn = x[k][j+1][i];
220: an = 0.5*(t0 + tn);
221: dn = PetscPowScalar(an,beta);
222: fn = dn*(tn - t0);
223: } else {
224: fn = zero;
225: }
227: if (k > 0) {
228: td = x[k-1][j][i];
229: ad = 0.5*(t0 + td);
230: dd = PetscPowScalar(ad,beta);
231: fd = dd*(t0 - td);
232: } else {
233: fd = zero;
234: }
236: if (k < mz-1) {
237: tu = x[k+1][j][i];
238: au = 0.5*(t0 + tu);
239: du = PetscPowScalar(au,beta);
240: fu = du*(tu - t0);
241: } else {
242: fu = zero;
243: }
245: } else if (i == mx-1) {
247: /* right-hand (east) boundary */
248: tw = x[k][j][i-1];
249: aw = 0.5*(t0 + tw);
250: dw = PetscPowScalar(aw,beta);
251: fw = dw*(t0 - tw);
253: te = tright;
254: ae = 0.5*(t0 + te);
255: de = PetscPowScalar(ae,beta);
256: fe = de*(te - t0);
258: if (j > 0) {
259: ts = x[k][j-1][i];
260: as = 0.5*(t0 + ts);
261: ds = PetscPowScalar(as,beta);
262: fs = ds*(t0 - ts);
263: } else {
264: fs = zero;
265: }
267: if (j < my-1) {
268: tn = x[k][j+1][i];
269: an = 0.5*(t0 + tn);
270: dn = PetscPowScalar(an,beta);
271: fn = dn*(tn - t0);
272: } else {
273: fn = zero;
274: }
276: if (k > 0) {
277: td = x[k-1][j][i];
278: ad = 0.5*(t0 + td);
279: dd = PetscPowScalar(ad,beta);
280: fd = dd*(t0 - td);
281: } else {
282: fd = zero;
283: }
285: if (k < mz-1) {
286: tu = x[k+1][j][i];
287: au = 0.5*(t0 + tu);
288: du = PetscPowScalar(au,beta);
289: fu = du*(tu - t0);
290: } else {
291: fu = zero;
292: }
294: } else if (j == 0) {
296: /* bottom (south) boundary, and i <> 0 or mx-1 */
297: tw = x[k][j][i-1];
298: aw = 0.5*(t0 + tw);
299: dw = PetscPowScalar(aw,beta);
300: fw = dw*(t0 - tw);
302: te = x[k][j][i+1];
303: ae = 0.5*(t0 + te);
304: de = PetscPowScalar(ae,beta);
305: fe = de*(te - t0);
307: fs = zero;
309: tn = x[k][j+1][i];
310: an = 0.5*(t0 + tn);
311: dn = PetscPowScalar(an,beta);
312: fn = dn*(tn - t0);
314: if (k > 0) {
315: td = x[k-1][j][i];
316: ad = 0.5*(t0 + td);
317: dd = PetscPowScalar(ad,beta);
318: fd = dd*(t0 - td);
319: } else {
320: fd = zero;
321: }
323: if (k < mz-1) {
324: tu = x[k+1][j][i];
325: au = 0.5*(t0 + tu);
326: du = PetscPowScalar(au,beta);
327: fu = du*(tu - t0);
328: } else {
329: fu = zero;
330: }
332: } else if (j == my-1) {
334: /* top (north) boundary, and i <> 0 or mx-1 */
335: tw = x[k][j][i-1];
336: aw = 0.5*(t0 + tw);
337: dw = PetscPowScalar(aw,beta);
338: fw = dw*(t0 - tw);
340: te = x[k][j][i+1];
341: ae = 0.5*(t0 + te);
342: de = PetscPowScalar(ae,beta);
343: fe = de*(te - t0);
345: ts = x[k][j-1][i];
346: as = 0.5*(t0 + ts);
347: ds = PetscPowScalar(as,beta);
348: fs = ds*(t0 - ts);
350: fn = zero;
352: if (k > 0) {
353: td = x[k-1][j][i];
354: ad = 0.5*(t0 + td);
355: dd = PetscPowScalar(ad,beta);
356: fd = dd*(t0 - td);
357: } else {
358: fd = zero;
359: }
361: if (k < mz-1) {
362: tu = x[k+1][j][i];
363: au = 0.5*(t0 + tu);
364: du = PetscPowScalar(au,beta);
365: fu = du*(tu - t0);
366: } else {
367: fu = zero;
368: }
370: } else if (k == 0) {
372: /* down boundary (interior only) */
373: tw = x[k][j][i-1];
374: aw = 0.5*(t0 + tw);
375: dw = PetscPowScalar(aw,beta);
376: fw = dw*(t0 - tw);
378: te = x[k][j][i+1];
379: ae = 0.5*(t0 + te);
380: de = PetscPowScalar(ae,beta);
381: fe = de*(te - t0);
383: ts = x[k][j-1][i];
384: as = 0.5*(t0 + ts);
385: ds = PetscPowScalar(as,beta);
386: fs = ds*(t0 - ts);
388: tn = x[k][j+1][i];
389: an = 0.5*(t0 + tn);
390: dn = PetscPowScalar(an,beta);
391: fn = dn*(tn - t0);
393: fd = zero;
395: tu = x[k+1][j][i];
396: au = 0.5*(t0 + tu);
397: du = PetscPowScalar(au,beta);
398: fu = du*(tu - t0);
400: } else if (k == mz-1) {
402: /* up boundary (interior only) */
403: tw = x[k][j][i-1];
404: aw = 0.5*(t0 + tw);
405: dw = PetscPowScalar(aw,beta);
406: fw = dw*(t0 - tw);
408: te = x[k][j][i+1];
409: ae = 0.5*(t0 + te);
410: de = PetscPowScalar(ae,beta);
411: fe = de*(te - t0);
413: ts = x[k][j-1][i];
414: as = 0.5*(t0 + ts);
415: ds = PetscPowScalar(as,beta);
416: fs = ds*(t0 - ts);
418: tn = x[k][j+1][i];
419: an = 0.5*(t0 + tn);
420: dn = PetscPowScalar(an,beta);
421: fn = dn*(tn - t0);
423: td = x[k-1][j][i];
424: ad = 0.5*(t0 + td);
425: dd = PetscPowScalar(ad,beta);
426: fd = dd*(t0 - td);
428: fu = zero;
429: }
431: f[k][j][i] = -hyhzdhx*(fe-fw) - hzhxdhy*(fn-fs) - hxhydhz*(fu-fd);
432: }
433: }
434: }
435: DMDAVecRestoreArray(da,localX,&x);
436: DMDAVecRestoreArray(da,F,&f);
437: DMRestoreLocalVector(da,&localX);
438: PetscLogFlops((22.0 + 4.0*POWFLOP)*ym*xm);
439: return 0;
440: }
441: /* -------------------- Evaluate Jacobian F(x) --------------------- */
442: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat jac,void *ptr)
443: {
444: AppCtx *user = (AppCtx*)ptr;
445: PetscInt i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm;
446: PetscScalar one = 1.0;
447: PetscScalar hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy;
448: PetscScalar t0,tn,ts,te,tw,an,as,ae,aw,dn,ds,de,dw;
449: PetscScalar tleft,tright,beta,td,ad,dd,tu,au,du,v[7],bm1,coef;
450: PetscScalar ***x,bn,bs,be,bw,bu,bd,gn,gs,ge,gw,gu,gd;
451: Vec localX;
452: MatStencil c[7],row;
453: DM da;
456: SNESGetDM(snes,&da);
457: DMGetLocalVector(da,&localX);
458: DMDAGetInfo(da,NULL,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
459: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1); hz = one/(PetscReal)(mz-1);
460: hxhydhz = hx*hy/hz; hyhzdhx = hy*hz/hx; hzhxdhy = hz*hx/hy;
461: tleft = user->tleft; tright = user->tright;
462: beta = user->beta; bm1 = user->bm1; coef = user->coef;
464: /* Get ghost points */
465: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
466: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
467: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
468: DMDAVecGetArray(da,localX,&x);
470: /* Evaluate Jacobian of function */
471: for (k=zs; k<zs+zm; k++) {
472: for (j=ys; j<ys+ym; j++) {
473: for (i=xs; i<xs+xm; i++) {
474: t0 = x[k][j][i];
475: row.k = k; row.j = j; row.i = i;
476: if (i > 0 && i < mx-1 && j > 0 && j < my-1 && k > 0 && k < mz-1) {
478: /* general interior volume */
480: tw = x[k][j][i-1];
481: aw = 0.5*(t0 + tw);
482: bw = PetscPowScalar(aw,bm1);
483: /* dw = bw * aw */
484: dw = PetscPowScalar(aw,beta);
485: gw = coef*bw*(t0 - tw);
487: te = x[k][j][i+1];
488: ae = 0.5*(t0 + te);
489: be = PetscPowScalar(ae,bm1);
490: /* de = be * ae; */
491: de = PetscPowScalar(ae,beta);
492: ge = coef*be*(te - t0);
494: ts = x[k][j-1][i];
495: as = 0.5*(t0 + ts);
496: bs = PetscPowScalar(as,bm1);
497: /* ds = bs * as; */
498: ds = PetscPowScalar(as,beta);
499: gs = coef*bs*(t0 - ts);
501: tn = x[k][j+1][i];
502: an = 0.5*(t0 + tn);
503: bn = PetscPowScalar(an,bm1);
504: /* dn = bn * an; */
505: dn = PetscPowScalar(an,beta);
506: gn = coef*bn*(tn - t0);
508: td = x[k-1][j][i];
509: ad = 0.5*(t0 + td);
510: bd = PetscPowScalar(ad,bm1);
511: /* dd = bd * ad; */
512: dd = PetscPowScalar(ad,beta);
513: gd = coef*bd*(t0 - td);
515: tu = x[k+1][j][i];
516: au = 0.5*(t0 + tu);
517: bu = PetscPowScalar(au,bm1);
518: /* du = bu * au; */
519: du = PetscPowScalar(au,beta);
520: gu = coef*bu*(tu - t0);
522: c[0].k = k-1; c[0].j = j; c[0].i = i; v[0] = -hxhydhz*(dd - gd);
523: c[1].k = k; c[1].j = j-1; c[1].i = i;
524: v[1] = -hzhxdhy*(ds - gs);
525: c[2].k = k; c[2].j = j; c[2].i = i-1;
526: v[2] = -hyhzdhx*(dw - gw);
527: c[3].k = k; c[3].j = j; c[3].i = i;
528: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
529: c[4].k = k; c[4].j = j; c[4].i = i+1;
530: v[4] = -hyhzdhx*(de + ge);
531: c[5].k = k; c[5].j = j+1; c[5].i = i;
532: v[5] = -hzhxdhy*(dn + gn);
533: c[6].k = k+1; c[6].j = j; c[6].i = i;
534: v[6] = -hxhydhz*(du + gu);
535: MatSetValuesStencil(jac,1,&row,7,c,v,INSERT_VALUES);
537: } else if (i == 0) {
539: /* left-hand plane boundary */
540: tw = tleft;
541: aw = 0.5*(t0 + tw);
542: bw = PetscPowScalar(aw,bm1);
543: /* dw = bw * aw */
544: dw = PetscPowScalar(aw,beta);
545: gw = coef*bw*(t0 - tw);
547: te = x[k][j][i+1];
548: ae = 0.5*(t0 + te);
549: be = PetscPowScalar(ae,bm1);
550: /* de = be * ae; */
551: de = PetscPowScalar(ae,beta);
552: ge = coef*be*(te - t0);
554: /* left-hand bottom edge */
555: if (j == 0) {
557: tn = x[k][j+1][i];
558: an = 0.5*(t0 + tn);
559: bn = PetscPowScalar(an,bm1);
560: /* dn = bn * an; */
561: dn = PetscPowScalar(an,beta);
562: gn = coef*bn*(tn - t0);
564: /* left-hand bottom down corner */
565: if (k == 0) {
567: tu = x[k+1][j][i];
568: au = 0.5*(t0 + tu);
569: bu = PetscPowScalar(au,bm1);
570: /* du = bu * au; */
571: du = PetscPowScalar(au,beta);
572: gu = coef*bu*(tu - t0);
574: c[0].k = k; c[0].j = j; c[0].i = i;
575: v[0] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
576: c[1].k = k; c[1].j = j; c[1].i = i+1;
577: v[1] = -hyhzdhx*(de + ge);
578: c[2].k = k; c[2].j = j+1; c[2].i = i;
579: v[2] = -hzhxdhy*(dn + gn);
580: c[3].k = k+1; c[3].j = j; c[3].i = i;
581: v[3] = -hxhydhz*(du + gu);
582: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
584: /* left-hand bottom interior edge */
585: } else if (k < mz-1) {
587: tu = x[k+1][j][i];
588: au = 0.5*(t0 + tu);
589: bu = PetscPowScalar(au,bm1);
590: /* du = bu * au; */
591: du = PetscPowScalar(au,beta);
592: gu = coef*bu*(tu - t0);
594: td = x[k-1][j][i];
595: ad = 0.5*(t0 + td);
596: bd = PetscPowScalar(ad,bm1);
597: /* dd = bd * ad; */
598: dd = PetscPowScalar(ad,beta);
599: gd = coef*bd*(td - t0);
601: c[0].k = k-1; c[0].j = j; c[0].i = i;
602: v[0] = -hxhydhz*(dd - gd);
603: c[1].k = k; c[1].j = j; c[1].i = i;
604: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
605: c[2].k = k; c[2].j = j; c[2].i = i+1;
606: v[2] = -hyhzdhx*(de + ge);
607: c[3].k = k; c[3].j = j+1; c[3].i = i;
608: v[3] = -hzhxdhy*(dn + gn);
609: c[4].k = k+1; c[4].j = j; c[4].i = i;
610: v[4] = -hxhydhz*(du + gu);
611: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
613: /* left-hand bottom up corner */
614: } else {
616: td = x[k-1][j][i];
617: ad = 0.5*(t0 + td);
618: bd = PetscPowScalar(ad,bm1);
619: /* dd = bd * ad; */
620: dd = PetscPowScalar(ad,beta);
621: gd = coef*bd*(td - t0);
623: c[0].k = k-1; c[0].j = j; c[0].i = i;
624: v[0] = -hxhydhz*(dd - gd);
625: c[1].k = k; c[1].j = j; c[1].i = i;
626: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
627: c[2].k = k; c[2].j = j; c[2].i = i+1;
628: v[2] = -hyhzdhx*(de + ge);
629: c[3].k = k; c[3].j = j+1; c[3].i = i;
630: v[3] = -hzhxdhy*(dn + gn);
631: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
632: }
634: /* left-hand top edge */
635: } else if (j == my-1) {
637: ts = x[k][j-1][i];
638: as = 0.5*(t0 + ts);
639: bs = PetscPowScalar(as,bm1);
640: /* ds = bs * as; */
641: ds = PetscPowScalar(as,beta);
642: gs = coef*bs*(ts - t0);
644: /* left-hand top down corner */
645: if (k == 0) {
647: tu = x[k+1][j][i];
648: au = 0.5*(t0 + tu);
649: bu = PetscPowScalar(au,bm1);
650: /* du = bu * au; */
651: du = PetscPowScalar(au,beta);
652: gu = coef*bu*(tu - t0);
654: c[0].k = k; c[0].j = j-1; c[0].i = i;
655: v[0] = -hzhxdhy*(ds - gs);
656: c[1].k = k; c[1].j = j; c[1].i = i;
657: v[1] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
658: c[2].k = k; c[2].j = j; c[2].i = i+1;
659: v[2] = -hyhzdhx*(de + ge);
660: c[3].k = k+1; c[3].j = j; c[3].i = i;
661: v[3] = -hxhydhz*(du + gu);
662: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
664: /* left-hand top interior edge */
665: } else if (k < mz-1) {
667: tu = x[k+1][j][i];
668: au = 0.5*(t0 + tu);
669: bu = PetscPowScalar(au,bm1);
670: /* du = bu * au; */
671: du = PetscPowScalar(au,beta);
672: gu = coef*bu*(tu - t0);
674: td = x[k-1][j][i];
675: ad = 0.5*(t0 + td);
676: bd = PetscPowScalar(ad,bm1);
677: /* dd = bd * ad; */
678: dd = PetscPowScalar(ad,beta);
679: gd = coef*bd*(td - t0);
681: c[0].k = k-1; c[0].j = j; c[0].i = i;
682: v[0] = -hxhydhz*(dd - gd);
683: c[1].k = k; c[1].j = j-1; c[1].i = i;
684: v[1] = -hzhxdhy*(ds - gs);
685: c[2].k = k; c[2].j = j; c[2].i = i;
686: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
687: c[3].k = k; c[3].j = j; c[3].i = i+1;
688: v[3] = -hyhzdhx*(de + ge);
689: c[4].k = k+1; c[4].j = j; c[4].i = i;
690: v[4] = -hxhydhz*(du + gu);
691: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
693: /* left-hand top up corner */
694: } else {
696: td = x[k-1][j][i];
697: ad = 0.5*(t0 + td);
698: bd = PetscPowScalar(ad,bm1);
699: /* dd = bd * ad; */
700: dd = PetscPowScalar(ad,beta);
701: gd = coef*bd*(td - t0);
703: c[0].k = k-1; c[0].j = j; c[0].i = i;
704: v[0] = -hxhydhz*(dd - gd);
705: c[1].k = k; c[1].j = j-1; c[1].i = i;
706: v[1] = -hzhxdhy*(ds - gs);
707: c[2].k = k; c[2].j = j; c[2].i = i;
708: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
709: c[3].k = k; c[3].j = j; c[3].i = i+1;
710: v[3] = -hyhzdhx*(de + ge);
711: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
712: }
714: } else {
716: ts = x[k][j-1][i];
717: as = 0.5*(t0 + ts);
718: bs = PetscPowScalar(as,bm1);
719: /* ds = bs * as; */
720: ds = PetscPowScalar(as,beta);
721: gs = coef*bs*(t0 - ts);
723: tn = x[k][j+1][i];
724: an = 0.5*(t0 + tn);
725: bn = PetscPowScalar(an,bm1);
726: /* dn = bn * an; */
727: dn = PetscPowScalar(an,beta);
728: gn = coef*bn*(tn - t0);
730: /* left-hand down interior edge */
731: if (k == 0) {
733: tu = x[k+1][j][i];
734: au = 0.5*(t0 + tu);
735: bu = PetscPowScalar(au,bm1);
736: /* du = bu * au; */
737: du = PetscPowScalar(au,beta);
738: gu = coef*bu*(tu - t0);
740: c[0].k = k; c[0].j = j-1; c[0].i = i;
741: v[0] = -hzhxdhy*(ds - gs);
742: c[1].k = k; c[1].j = j; c[1].i = i;
743: v[1] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
744: c[2].k = k; c[2].j = j; c[2].i = i+1;
745: v[2] = -hyhzdhx*(de + ge);
746: c[3].k = k; c[3].j = j+1; c[3].i = i;
747: v[3] = -hzhxdhy*(dn + gn);
748: c[4].k = k+1; c[4].j = j; c[4].i = i;
749: v[4] = -hxhydhz*(du + gu);
750: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
752: } else if (k == mz-1) { /* left-hand up interior edge */
754: td = x[k-1][j][i];
755: ad = 0.5*(t0 + td);
756: bd = PetscPowScalar(ad,bm1);
757: /* dd = bd * ad; */
758: dd = PetscPowScalar(ad,beta);
759: gd = coef*bd*(t0 - td);
761: c[0].k = k-1; c[0].j = j; c[0].i = i;
762: v[0] = -hxhydhz*(dd - gd);
763: c[1].k = k; c[1].j = j-1; c[1].i = i;
764: v[1] = -hzhxdhy*(ds - gs);
765: c[2].k = k; c[2].j = j; c[2].i = i;
766: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
767: c[3].k = k; c[3].j = j; c[3].i = i+1;
768: v[3] = -hyhzdhx*(de + ge);
769: c[4].k = k; c[4].j = j+1; c[4].i = i;
770: v[4] = -hzhxdhy*(dn + gn);
771: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
772: } else { /* left-hand interior plane */
774: td = x[k-1][j][i];
775: ad = 0.5*(t0 + td);
776: bd = PetscPowScalar(ad,bm1);
777: /* dd = bd * ad; */
778: dd = PetscPowScalar(ad,beta);
779: gd = coef*bd*(t0 - td);
781: tu = x[k+1][j][i];
782: au = 0.5*(t0 + tu);
783: bu = PetscPowScalar(au,bm1);
784: /* du = bu * au; */
785: du = PetscPowScalar(au,beta);
786: gu = coef*bu*(tu - t0);
788: c[0].k = k-1; c[0].j = j; c[0].i = i;
789: v[0] = -hxhydhz*(dd - gd);
790: c[1].k = k; c[1].j = j-1; c[1].i = i;
791: v[1] = -hzhxdhy*(ds - gs);
792: c[2].k = k; c[2].j = j; c[2].i = i;
793: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
794: c[3].k = k; c[3].j = j; c[3].i = i+1;
795: v[3] = -hyhzdhx*(de + ge);
796: c[4].k = k; c[4].j = j+1; c[4].i = i;
797: v[4] = -hzhxdhy*(dn + gn);
798: c[5].k = k+1; c[5].j = j; c[5].i = i;
799: v[5] = -hxhydhz*(du + gu);
800: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
801: }
802: }
804: } else if (i == mx-1) {
806: /* right-hand plane boundary */
807: tw = x[k][j][i-1];
808: aw = 0.5*(t0 + tw);
809: bw = PetscPowScalar(aw,bm1);
810: /* dw = bw * aw */
811: dw = PetscPowScalar(aw,beta);
812: gw = coef*bw*(t0 - tw);
814: te = tright;
815: ae = 0.5*(t0 + te);
816: be = PetscPowScalar(ae,bm1);
817: /* de = be * ae; */
818: de = PetscPowScalar(ae,beta);
819: ge = coef*be*(te - t0);
821: /* right-hand bottom edge */
822: if (j == 0) {
824: tn = x[k][j+1][i];
825: an = 0.5*(t0 + tn);
826: bn = PetscPowScalar(an,bm1);
827: /* dn = bn * an; */
828: dn = PetscPowScalar(an,beta);
829: gn = coef*bn*(tn - t0);
831: /* right-hand bottom down corner */
832: if (k == 0) {
834: tu = x[k+1][j][i];
835: au = 0.5*(t0 + tu);
836: bu = PetscPowScalar(au,bm1);
837: /* du = bu * au; */
838: du = PetscPowScalar(au,beta);
839: gu = coef*bu*(tu - t0);
841: c[0].k = k; c[0].j = j; c[0].i = i-1;
842: v[0] = -hyhzdhx*(dw - gw);
843: c[1].k = k; c[1].j = j; c[1].i = i;
844: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
845: c[2].k = k; c[2].j = j+1; c[2].i = i;
846: v[2] = -hzhxdhy*(dn + gn);
847: c[3].k = k+1; c[3].j = j; c[3].i = i;
848: v[3] = -hxhydhz*(du + gu);
849: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
851: /* right-hand bottom interior edge */
852: } else if (k < mz-1) {
854: tu = x[k+1][j][i];
855: au = 0.5*(t0 + tu);
856: bu = PetscPowScalar(au,bm1);
857: /* du = bu * au; */
858: du = PetscPowScalar(au,beta);
859: gu = coef*bu*(tu - t0);
861: td = x[k-1][j][i];
862: ad = 0.5*(t0 + td);
863: bd = PetscPowScalar(ad,bm1);
864: /* dd = bd * ad; */
865: dd = PetscPowScalar(ad,beta);
866: gd = coef*bd*(td - t0);
868: c[0].k = k-1; c[0].j = j; c[0].i = i;
869: v[0] = -hxhydhz*(dd - gd);
870: c[1].k = k; c[1].j = j; c[1].i = i-1;
871: v[1] = -hyhzdhx*(dw - gw);
872: c[2].k = k; c[2].j = j; c[2].i = i;
873: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
874: c[3].k = k; c[3].j = j+1; c[3].i = i;
875: v[3] = -hzhxdhy*(dn + gn);
876: c[4].k = k+1; c[4].j = j; c[4].i = i;
877: v[4] = -hxhydhz*(du + gu);
878: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
880: /* right-hand bottom up corner */
881: } else {
883: td = x[k-1][j][i];
884: ad = 0.5*(t0 + td);
885: bd = PetscPowScalar(ad,bm1);
886: /* dd = bd * ad; */
887: dd = PetscPowScalar(ad,beta);
888: gd = coef*bd*(td - t0);
890: c[0].k = k-1; c[0].j = j; c[0].i = i;
891: v[0] = -hxhydhz*(dd - gd);
892: c[1].k = k; c[1].j = j; c[1].i = i-1;
893: v[1] = -hyhzdhx*(dw - gw);
894: c[2].k = k; c[2].j = j; c[2].i = i;
895: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
896: c[3].k = k; c[3].j = j+1; c[3].i = i;
897: v[3] = -hzhxdhy*(dn + gn);
898: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
899: }
901: /* right-hand top edge */
902: } else if (j == my-1) {
904: ts = x[k][j-1][i];
905: as = 0.5*(t0 + ts);
906: bs = PetscPowScalar(as,bm1);
907: /* ds = bs * as; */
908: ds = PetscPowScalar(as,beta);
909: gs = coef*bs*(ts - t0);
911: /* right-hand top down corner */
912: if (k == 0) {
914: tu = x[k+1][j][i];
915: au = 0.5*(t0 + tu);
916: bu = PetscPowScalar(au,bm1);
917: /* du = bu * au; */
918: du = PetscPowScalar(au,beta);
919: gu = coef*bu*(tu - t0);
921: c[0].k = k; c[0].j = j-1; c[0].i = i;
922: v[0] = -hzhxdhy*(ds - gs);
923: c[1].k = k; c[1].j = j; c[1].i = i-1;
924: v[1] = -hyhzdhx*(dw - gw);
925: c[2].k = k; c[2].j = j; c[2].i = i;
926: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
927: c[3].k = k+1; c[3].j = j; c[3].i = i;
928: v[3] = -hxhydhz*(du + gu);
929: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
931: /* right-hand top interior edge */
932: } else if (k < mz-1) {
934: tu = x[k+1][j][i];
935: au = 0.5*(t0 + tu);
936: bu = PetscPowScalar(au,bm1);
937: /* du = bu * au; */
938: du = PetscPowScalar(au,beta);
939: gu = coef*bu*(tu - t0);
941: td = x[k-1][j][i];
942: ad = 0.5*(t0 + td);
943: bd = PetscPowScalar(ad,bm1);
944: /* dd = bd * ad; */
945: dd = PetscPowScalar(ad,beta);
946: gd = coef*bd*(td - t0);
948: c[0].k = k-1; c[0].j = j; c[0].i = i;
949: v[0] = -hxhydhz*(dd - gd);
950: c[1].k = k; c[1].j = j-1; c[1].i = i;
951: v[1] = -hzhxdhy*(ds - gs);
952: c[2].k = k; c[2].j = j; c[2].i = i-1;
953: v[2] = -hyhzdhx*(dw - gw);
954: c[3].k = k; c[3].j = j; c[3].i = i;
955: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
956: c[4].k = k+1; c[4].j = j; c[4].i = i;
957: v[4] = -hxhydhz*(du + gu);
958: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
960: /* right-hand top up corner */
961: } else {
963: td = x[k-1][j][i];
964: ad = 0.5*(t0 + td);
965: bd = PetscPowScalar(ad,bm1);
966: /* dd = bd * ad; */
967: dd = PetscPowScalar(ad,beta);
968: gd = coef*bd*(td - t0);
970: c[0].k = k-1; c[0].j = j; c[0].i = i;
971: v[0] = -hxhydhz*(dd - gd);
972: c[1].k = k; c[1].j = j-1; c[1].i = i;
973: v[1] = -hzhxdhy*(ds - gs);
974: c[2].k = k; c[2].j = j; c[2].i = i-1;
975: v[2] = -hyhzdhx*(dw - gw);
976: c[3].k = k; c[3].j = j; c[3].i = i;
977: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
978: MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES);
979: }
981: } else {
983: ts = x[k][j-1][i];
984: as = 0.5*(t0 + ts);
985: bs = PetscPowScalar(as,bm1);
986: /* ds = bs * as; */
987: ds = PetscPowScalar(as,beta);
988: gs = coef*bs*(t0 - ts);
990: tn = x[k][j+1][i];
991: an = 0.5*(t0 + tn);
992: bn = PetscPowScalar(an,bm1);
993: /* dn = bn * an; */
994: dn = PetscPowScalar(an,beta);
995: gn = coef*bn*(tn - t0);
997: /* right-hand down interior edge */
998: if (k == 0) {
1000: tu = x[k+1][j][i];
1001: au = 0.5*(t0 + tu);
1002: bu = PetscPowScalar(au,bm1);
1003: /* du = bu * au; */
1004: du = PetscPowScalar(au,beta);
1005: gu = coef*bu*(tu - t0);
1007: c[0].k = k; c[0].j = j-1; c[0].i = i;
1008: v[0] = -hzhxdhy*(ds - gs);
1009: c[1].k = k; c[1].j = j; c[1].i = i-1;
1010: v[1] = -hyhzdhx*(dw - gw);
1011: c[2].k = k; c[2].j = j; c[2].i = i;
1012: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1013: c[3].k = k; c[3].j = j+1; c[3].i = i;
1014: v[3] = -hzhxdhy*(dn + gn);
1015: c[4].k = k+1; c[4].j = j; c[4].i = i;
1016: v[4] = -hxhydhz*(du + gu);
1017: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1019: } else if (k == mz-1) { /* right-hand up interior edge */
1021: td = x[k-1][j][i];
1022: ad = 0.5*(t0 + td);
1023: bd = PetscPowScalar(ad,bm1);
1024: /* dd = bd * ad; */
1025: dd = PetscPowScalar(ad,beta);
1026: gd = coef*bd*(t0 - td);
1028: c[0].k = k-1; c[0].j = j; c[0].i = i;
1029: v[0] = -hxhydhz*(dd - gd);
1030: c[1].k = k; c[1].j = j-1; c[1].i = i;
1031: v[1] = -hzhxdhy*(ds - gs);
1032: c[2].k = k; c[2].j = j; c[2].i = i-1;
1033: v[2] = -hyhzdhx*(dw - gw);
1034: c[3].k = k; c[3].j = j; c[3].i = i;
1035: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1036: c[4].k = k; c[4].j = j+1; c[4].i = i;
1037: v[4] = -hzhxdhy*(dn + gn);
1038: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1040: } else { /* right-hand interior plane */
1042: td = x[k-1][j][i];
1043: ad = 0.5*(t0 + td);
1044: bd = PetscPowScalar(ad,bm1);
1045: /* dd = bd * ad; */
1046: dd = PetscPowScalar(ad,beta);
1047: gd = coef*bd*(t0 - td);
1049: tu = x[k+1][j][i];
1050: au = 0.5*(t0 + tu);
1051: bu = PetscPowScalar(au,bm1);
1052: /* du = bu * au; */
1053: du = PetscPowScalar(au,beta);
1054: gu = coef*bu*(tu - t0);
1056: c[0].k = k-1; c[0].j = j; c[0].i = i;
1057: v[0] = -hxhydhz*(dd - gd);
1058: c[1].k = k; c[1].j = j-1; c[1].i = i;
1059: v[1] = -hzhxdhy*(ds - gs);
1060: c[2].k = k; c[2].j = j; c[2].i = i-1;
1061: v[2] = -hyhzdhx*(dw - gw);
1062: c[3].k = k; c[3].j = j; c[3].i = i;
1063: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1064: c[4].k = k; c[4].j = j+1; c[4].i = i;
1065: v[4] = -hzhxdhy*(dn + gn);
1066: c[5].k = k+1; c[5].j = j; c[5].i = i;
1067: v[5] = -hxhydhz*(du + gu);
1068: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1069: }
1070: }
1072: } else if (j == 0) {
1074: tw = x[k][j][i-1];
1075: aw = 0.5*(t0 + tw);
1076: bw = PetscPowScalar(aw,bm1);
1077: /* dw = bw * aw */
1078: dw = PetscPowScalar(aw,beta);
1079: gw = coef*bw*(t0 - tw);
1081: te = x[k][j][i+1];
1082: ae = 0.5*(t0 + te);
1083: be = PetscPowScalar(ae,bm1);
1084: /* de = be * ae; */
1085: de = PetscPowScalar(ae,beta);
1086: ge = coef*be*(te - t0);
1088: tn = x[k][j+1][i];
1089: an = 0.5*(t0 + tn);
1090: bn = PetscPowScalar(an,bm1);
1091: /* dn = bn * an; */
1092: dn = PetscPowScalar(an,beta);
1093: gn = coef*bn*(tn - t0);
1095: /* bottom down interior edge */
1096: if (k == 0) {
1098: tu = x[k+1][j][i];
1099: au = 0.5*(t0 + tu);
1100: bu = PetscPowScalar(au,bm1);
1101: /* du = bu * au; */
1102: du = PetscPowScalar(au,beta);
1103: gu = coef*bu*(tu - t0);
1105: c[0].k = k; c[0].j = j; c[0].i = i-1;
1106: v[0] = -hyhzdhx*(dw - gw);
1107: c[1].k = k; c[1].j = j; c[1].i = i;
1108: v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1109: c[2].k = k; c[2].j = j; c[2].i = i+1;
1110: v[2] = -hyhzdhx*(de + ge);
1111: c[3].k = k; c[3].j = j+1; c[3].i = i;
1112: v[3] = -hzhxdhy*(dn + gn);
1113: c[4].k = k+1; c[4].j = j; c[4].i = i;
1114: v[4] = -hxhydhz*(du + gu);
1115: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1117: } else if (k == mz-1) { /* bottom up interior edge */
1119: td = x[k-1][j][i];
1120: ad = 0.5*(t0 + td);
1121: bd = PetscPowScalar(ad,bm1);
1122: /* dd = bd * ad; */
1123: dd = PetscPowScalar(ad,beta);
1124: gd = coef*bd*(td - t0);
1126: c[0].k = k-1; c[0].j = j; c[0].i = i;
1127: v[0] = -hxhydhz*(dd - gd);
1128: c[1].k = k; c[1].j = j; c[1].i = i-1;
1129: v[1] = -hyhzdhx*(dw - gw);
1130: c[2].k = k; c[2].j = j; c[2].i = i;
1131: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1132: c[3].k = k; c[3].j = j; c[3].i = i+1;
1133: v[3] = -hyhzdhx*(de + ge);
1134: c[4].k = k; c[4].j = j+1; c[4].i = i;
1135: v[4] = -hzhxdhy*(dn + gn);
1136: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1138: } else { /* bottom interior plane */
1140: tu = x[k+1][j][i];
1141: au = 0.5*(t0 + tu);
1142: bu = PetscPowScalar(au,bm1);
1143: /* du = bu * au; */
1144: du = PetscPowScalar(au,beta);
1145: gu = coef*bu*(tu - t0);
1147: td = x[k-1][j][i];
1148: ad = 0.5*(t0 + td);
1149: bd = PetscPowScalar(ad,bm1);
1150: /* dd = bd * ad; */
1151: dd = PetscPowScalar(ad,beta);
1152: gd = coef*bd*(td - t0);
1154: c[0].k = k-1; c[0].j = j; c[0].i = i;
1155: v[0] = -hxhydhz*(dd - gd);
1156: c[1].k = k; c[1].j = j; c[1].i = i-1;
1157: v[1] = -hyhzdhx*(dw - gw);
1158: c[2].k = k; c[2].j = j; c[2].i = i;
1159: v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1160: c[3].k = k; c[3].j = j; c[3].i = i+1;
1161: v[3] = -hyhzdhx*(de + ge);
1162: c[4].k = k; c[4].j = j+1; c[4].i = i;
1163: v[4] = -hzhxdhy*(dn + gn);
1164: c[5].k = k+1; c[5].j = j; c[5].i = i;
1165: v[5] = -hxhydhz*(du + gu);
1166: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1167: }
1169: } else if (j == my-1) {
1171: tw = x[k][j][i-1];
1172: aw = 0.5*(t0 + tw);
1173: bw = PetscPowScalar(aw,bm1);
1174: /* dw = bw * aw */
1175: dw = PetscPowScalar(aw,beta);
1176: gw = coef*bw*(t0 - tw);
1178: te = x[k][j][i+1];
1179: ae = 0.5*(t0 + te);
1180: be = PetscPowScalar(ae,bm1);
1181: /* de = be * ae; */
1182: de = PetscPowScalar(ae,beta);
1183: ge = coef*be*(te - t0);
1185: ts = x[k][j-1][i];
1186: as = 0.5*(t0 + ts);
1187: bs = PetscPowScalar(as,bm1);
1188: /* ds = bs * as; */
1189: ds = PetscPowScalar(as,beta);
1190: gs = coef*bs*(t0 - ts);
1192: /* top down interior edge */
1193: if (k == 0) {
1195: tu = x[k+1][j][i];
1196: au = 0.5*(t0 + tu);
1197: bu = PetscPowScalar(au,bm1);
1198: /* du = bu * au; */
1199: du = PetscPowScalar(au,beta);
1200: gu = coef*bu*(tu - t0);
1202: c[0].k = k; c[0].j = j-1; c[0].i = i;
1203: v[0] = -hzhxdhy*(ds - gs);
1204: c[1].k = k; c[1].j = j; c[1].i = i-1;
1205: v[1] = -hyhzdhx*(dw - gw);
1206: c[2].k = k; c[2].j = j; c[2].i = i;
1207: v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1208: c[3].k = k; c[3].j = j; c[3].i = i+1;
1209: v[3] = -hyhzdhx*(de + ge);
1210: c[4].k = k+1; c[4].j = j; c[4].i = i;
1211: v[4] = -hxhydhz*(du + gu);
1212: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1214: } else if (k == mz-1) { /* top up interior edge */
1216: td = x[k-1][j][i];
1217: ad = 0.5*(t0 + td);
1218: bd = PetscPowScalar(ad,bm1);
1219: /* dd = bd * ad; */
1220: dd = PetscPowScalar(ad,beta);
1221: gd = coef*bd*(td - t0);
1223: c[0].k = k-1; c[0].j = j; c[0].i = i;
1224: v[0] = -hxhydhz*(dd - gd);
1225: c[1].k = k; c[1].j = j-1; c[1].i = i;
1226: v[1] = -hzhxdhy*(ds - gs);
1227: c[2].k = k; c[2].j = j; c[2].i = i-1;
1228: v[2] = -hyhzdhx*(dw - gw);
1229: c[3].k = k; c[3].j = j; c[3].i = i;
1230: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1231: c[4].k = k; c[4].j = j; c[4].i = i+1;
1232: v[4] = -hyhzdhx*(de + ge);
1233: MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES);
1235: } else { /* top interior plane */
1237: tu = x[k+1][j][i];
1238: au = 0.5*(t0 + tu);
1239: bu = PetscPowScalar(au,bm1);
1240: /* du = bu * au; */
1241: du = PetscPowScalar(au,beta);
1242: gu = coef*bu*(tu - t0);
1244: td = x[k-1][j][i];
1245: ad = 0.5*(t0 + td);
1246: bd = PetscPowScalar(ad,bm1);
1247: /* dd = bd * ad; */
1248: dd = PetscPowScalar(ad,beta);
1249: gd = coef*bd*(td - t0);
1251: c[0].k = k-1; c[0].j = j; c[0].i = i;
1252: v[0] = -hxhydhz*(dd - gd);
1253: c[1].k = k; c[1].j = j-1; c[1].i = i;
1254: v[1] = -hzhxdhy*(ds - gs);
1255: c[2].k = k; c[2].j = j; c[2].i = i-1;
1256: v[2] = -hyhzdhx*(dw - gw);
1257: c[3].k = k; c[3].j = j; c[3].i = i;
1258: v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu);
1259: c[4].k = k; c[4].j = j; c[4].i = i+1;
1260: v[4] = -hyhzdhx*(de + ge);
1261: c[5].k = k+1; c[5].j = j; c[5].i = i;
1262: v[5] = -hxhydhz*(du + gu);
1263: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1264: }
1266: } else if (k == 0) {
1268: /* down interior plane */
1270: tw = x[k][j][i-1];
1271: aw = 0.5*(t0 + tw);
1272: bw = PetscPowScalar(aw,bm1);
1273: /* dw = bw * aw */
1274: dw = PetscPowScalar(aw,beta);
1275: gw = coef*bw*(t0 - tw);
1277: te = x[k][j][i+1];
1278: ae = 0.5*(t0 + te);
1279: be = PetscPowScalar(ae,bm1);
1280: /* de = be * ae; */
1281: de = PetscPowScalar(ae,beta);
1282: ge = coef*be*(te - t0);
1284: ts = x[k][j-1][i];
1285: as = 0.5*(t0 + ts);
1286: bs = PetscPowScalar(as,bm1);
1287: /* ds = bs * as; */
1288: ds = PetscPowScalar(as,beta);
1289: gs = coef*bs*(t0 - ts);
1291: tn = x[k][j+1][i];
1292: an = 0.5*(t0 + tn);
1293: bn = PetscPowScalar(an,bm1);
1294: /* dn = bn * an; */
1295: dn = PetscPowScalar(an,beta);
1296: gn = coef*bn*(tn - t0);
1298: tu = x[k+1][j][i];
1299: au = 0.5*(t0 + tu);
1300: bu = PetscPowScalar(au,bm1);
1301: /* du = bu * au; */
1302: du = PetscPowScalar(au,beta);
1303: gu = coef*bu*(tu - t0);
1305: c[0].k = k; c[0].j = j-1; c[0].i = i;
1306: v[0] = -hzhxdhy*(ds - gs);
1307: c[1].k = k; c[1].j = j; c[1].i = i-1;
1308: v[1] = -hyhzdhx*(dw - gw);
1309: c[2].k = k; c[2].j = j; c[2].i = i;
1310: v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu);
1311: c[3].k = k; c[3].j = j; c[3].i = i+1;
1312: v[3] = -hyhzdhx*(de + ge);
1313: c[4].k = k; c[4].j = j+1; c[4].i = i;
1314: v[4] = -hzhxdhy*(dn + gn);
1315: c[5].k = k+1; c[5].j = j; c[5].i = i;
1316: v[5] = -hxhydhz*(du + gu);
1317: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1319: } else if (k == mz-1) {
1321: /* up interior plane */
1323: tw = x[k][j][i-1];
1324: aw = 0.5*(t0 + tw);
1325: bw = PetscPowScalar(aw,bm1);
1326: /* dw = bw * aw */
1327: dw = PetscPowScalar(aw,beta);
1328: gw = coef*bw*(t0 - tw);
1330: te = x[k][j][i+1];
1331: ae = 0.5*(t0 + te);
1332: be = PetscPowScalar(ae,bm1);
1333: /* de = be * ae; */
1334: de = PetscPowScalar(ae,beta);
1335: ge = coef*be*(te - t0);
1337: ts = x[k][j-1][i];
1338: as = 0.5*(t0 + ts);
1339: bs = PetscPowScalar(as,bm1);
1340: /* ds = bs * as; */
1341: ds = PetscPowScalar(as,beta);
1342: gs = coef*bs*(t0 - ts);
1344: tn = x[k][j+1][i];
1345: an = 0.5*(t0 + tn);
1346: bn = PetscPowScalar(an,bm1);
1347: /* dn = bn * an; */
1348: dn = PetscPowScalar(an,beta);
1349: gn = coef*bn*(tn - t0);
1351: td = x[k-1][j][i];
1352: ad = 0.5*(t0 + td);
1353: bd = PetscPowScalar(ad,bm1);
1354: /* dd = bd * ad; */
1355: dd = PetscPowScalar(ad,beta);
1356: gd = coef*bd*(t0 - td);
1358: c[0].k = k-1; c[0].j = j; c[0].i = i;
1359: v[0] = -hxhydhz*(dd - gd);
1360: c[1].k = k; c[1].j = j-1; c[1].i = i;
1361: v[1] = -hzhxdhy*(ds - gs);
1362: c[2].k = k; c[2].j = j; c[2].i = i-1;
1363: v[2] = -hyhzdhx*(dw - gw);
1364: c[3].k = k; c[3].j = j; c[3].i = i;
1365: v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd);
1366: c[4].k = k; c[4].j = j; c[4].i = i+1;
1367: v[4] = -hyhzdhx*(de + ge);
1368: c[5].k = k; c[5].j = j+1; c[5].i = i;
1369: v[5] = -hzhxdhy*(dn + gn);
1370: MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES);
1371: }
1372: }
1373: }
1374: }
1375: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
1376: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
1377: if (jac != J) {
1378: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1379: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1380: }
1381: DMDAVecRestoreArray(da,localX,&x);
1382: DMRestoreLocalVector(da,&localX);
1384: PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym);
1385: return 0;
1386: }
1388: /*TEST
1390: test:
1391: nsize: 4
1392: args: -snes_monitor_short -pc_mg_type full -ksp_type fgmres -pc_type mg -snes_view -pc_mg_levels 2 -pc_mg_galerkin pmat
1393: requires: !single
1395: TEST*/