Actual source code: ex18.c
petsc-3.3-p7 2013-05-11
2: static char help[] ="Nonlinear Radiative Transport PDE with multigrid in 2d.\n\
3: Uses 2-dimensional distributed arrays.\n\
4: A 2-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: /*T
15: Concepts: SNES^solving a system of nonlinear equations
16: Concepts: DMDA^using distributed arrays
17: Concepts: multigrid;
18: Processors: n
19: T*/
21: /*
22:
23: This example models the partial differential equation
24:
25: - Div(alpha* T^beta (GRAD T)) = 0.
26:
27: where beta = 2.5 and alpha = 1.0
28:
29: BC: T_left = 1.0, T_right = 0.1, dT/dn_top = dTdn_bottom = 0.
30:
31: in the unit square, which is uniformly discretized in each of x and
32: y in this simple encoding. The degrees of freedom are cell centered.
33:
34: A finite volume approximation with the usual 5-point stencil
35: is used to discretize the boundary value problem to obtain a
36: nonlinear system of equations.
38: This code was contributed by David Keyes
39:
40: */
42: #include <petscsnes.h>
43: #include <petscdmda.h>
44: #include <petscpcmg.h>
46: /* User-defined application context */
48: typedef struct {
49: PetscReal tleft,tright; /* Dirichlet boundary conditions */
50: PetscReal beta,bm1,coef; /* nonlinear diffusivity parameterizations */
51: } AppCtx;
53: #define POWFLOP 5 /* assume a pow() takes five flops */
55: extern PetscErrorCode FormInitialGuess(DM,Vec);
56: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
57: extern PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
61: int main(int argc,char **argv)
62: {
63: SNES snes;
64: AppCtx user;
66: PetscInt its,lits;
67: PetscReal litspit;
68: DM da;
70: PetscInitialize(&argc,&argv,PETSC_NULL,help);
72: /* set problem parameters */
73: user.tleft = 1.0;
74: user.tright = 0.1;
75: user.beta = 2.5;
76: PetscOptionsGetReal(PETSC_NULL,"-tleft",&user.tleft,PETSC_NULL);
77: PetscOptionsGetReal(PETSC_NULL,"-tright",&user.tright,PETSC_NULL);
78: PetscOptionsGetReal(PETSC_NULL,"-beta",&user.beta,PETSC_NULL);
79: user.bm1 = user.beta - 1.0;
80: user.coef = user.beta/2.0;
83: /*
84: Create the multilevel DM data structure
85: */
86: SNESCreate(PETSC_COMM_WORLD,&snes);
88: /*
89: Set the DMDA (grid structure) for the grids.
90: */
91: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-5,-5,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);
92: DMSetApplicationContext(da,&user);
93: SNESSetDM(snes,(DM)da);
95: /*
96: Create the nonlinear solver, and tell it the functions to use
97: */
98: SNESSetFunction(snes,PETSC_NULL,FormFunction,&user);
99: SNESSetJacobian(snes,PETSC_NULL,PETSC_NULL,FormJacobian,&user);
100: SNESSetFromOptions(snes);
101: DMSetInitialGuess(da,FormInitialGuess);
103: SNESSolve(snes,PETSC_NULL,PETSC_NULL);
104: SNESGetIterationNumber(snes,&its);
105: SNESGetLinearSolveIterations(snes,&lits);
106: litspit = ((PetscReal)lits)/((PetscReal)its);
107: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
108: PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %D\n",lits);
109: PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / SNES = %e\n",litspit);
111: DMDestroy(&da);
112: SNESDestroy(&snes);
113: PetscFinalize();
115: return 0;
116: }
117: /* -------------------- Form initial approximation ----------------- */
120: PetscErrorCode FormInitialGuess(DM da,Vec X)
121: {
122: AppCtx *user;
123: PetscInt i,j,xs,ys,xm,ym;
125: PetscReal tleft;
126: PetscScalar **x;
129: DMGetApplicationContext(da,&user);
130: tleft = user->tleft;
131: /* Get ghost points */
132: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
133: DMDAVecGetArray(da,X,&x);
135: /* Compute initial guess */
136: for (j=ys; j<ys+ym; j++) {
137: for (i=xs; i<xs+xm; i++) {
138: x[j][i] = tleft;
139: }
140: }
141: DMDAVecRestoreArray(da,X,&x);
142: return(0);
143: }
144: /* -------------------- Evaluate Function F(x) --------------------- */
147: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void* ptr)
148: {
149: AppCtx *user = (AppCtx*)ptr;
151: PetscInt i,j,mx,my,xs,ys,xm,ym;
152: PetscScalar zero = 0.0,one = 1.0;
153: PetscScalar hx,hy,hxdhy,hydhx;
154: 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;
155: PetscScalar tleft,tright,beta;
156: PetscScalar **x,**f;
157: Vec localX;
158: DM da;
161: SNESGetDM(snes,&da);
162: DMGetLocalVector(da,&localX);
163: DMDAGetInfo(da,PETSC_NULL,&mx,&my,0,0,0,0,0,0,0,0,0,0);
164: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1);
165: hxdhy = hx/hy; hydhx = hy/hx;
166: tleft = user->tleft; tright = user->tright;
167: beta = user->beta;
168:
169: /* Get ghost points */
170: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
171: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
172: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
173: DMDAVecGetArray(da,localX,&x);
174: DMDAVecGetArray(da,F,&f);
176: /* Evaluate function */
177: for (j=ys; j<ys+ym; j++) {
178: for (i=xs; i<xs+xm; i++) {
179: t0 = x[j][i];
181: if (i > 0 && i < mx-1 && j > 0 && j < my-1) {
183: /* general interior volume */
185: tw = x[j][i-1];
186: aw = 0.5*(t0 + tw);
187: dw = PetscPowScalar(aw,beta);
188: fw = dw*(t0 - tw);
190: te = x[j][i+1];
191: ae = 0.5*(t0 + te);
192: de = PetscPowScalar(ae,beta);
193: fe = de*(te - t0);
195: ts = x[j-1][i];
196: as = 0.5*(t0 + ts);
197: ds = PetscPowScalar(as,beta);
198: fs = ds*(t0 - ts);
199:
200: tn = x[j+1][i];
201: an = 0.5*(t0 + tn);
202: dn = PetscPowScalar(an,beta);
203: fn = dn*(tn - t0);
205: } else if (i == 0) {
207: /* left-hand boundary */
208: tw = tleft;
209: aw = 0.5*(t0 + tw);
210: dw = PetscPowScalar(aw,beta);
211: fw = dw*(t0 - tw);
213: te = x[j][i+1];
214: ae = 0.5*(t0 + te);
215: de = PetscPowScalar(ae,beta);
216: fe = de*(te - t0);
218: if (j > 0) {
219: ts = x[j-1][i];
220: as = 0.5*(t0 + ts);
221: ds = PetscPowScalar(as,beta);
222: fs = ds*(t0 - ts);
223: } else {
224: fs = zero;
225: }
227: if (j < my-1) {
228: tn = x[j+1][i];
229: an = 0.5*(t0 + tn);
230: dn = PetscPowScalar(an,beta);
231: fn = dn*(tn - t0);
232: } else {
233: fn = zero;
234: }
236: } else if (i == mx-1) {
238: /* right-hand boundary */
239: tw = x[j][i-1];
240: aw = 0.5*(t0 + tw);
241: dw = PetscPowScalar(aw,beta);
242: fw = dw*(t0 - tw);
243:
244: te = tright;
245: ae = 0.5*(t0 + te);
246: de = PetscPowScalar(ae,beta);
247: fe = de*(te - t0);
248:
249: if (j > 0) {
250: ts = x[j-1][i];
251: as = 0.5*(t0 + ts);
252: ds = PetscPowScalar(as,beta);
253: fs = ds*(t0 - ts);
254: } else {
255: fs = zero;
256: }
257:
258: if (j < my-1) {
259: tn = x[j+1][i];
260: an = 0.5*(t0 + tn);
261: dn = PetscPowScalar(an,beta);
262: fn = dn*(tn - t0);
263: } else {
264: fn = zero;
265: }
267: } else if (j == 0) {
269: /* bottom boundary,and i <> 0 or mx-1 */
270: tw = x[j][i-1];
271: aw = 0.5*(t0 + tw);
272: dw = PetscPowScalar(aw,beta);
273: fw = dw*(t0 - tw);
275: te = x[j][i+1];
276: ae = 0.5*(t0 + te);
277: de = PetscPowScalar(ae,beta);
278: fe = de*(te - t0);
280: fs = zero;
282: tn = x[j+1][i];
283: an = 0.5*(t0 + tn);
284: dn = PetscPowScalar(an,beta);
285: fn = dn*(tn - t0);
287: } else if (j == my-1) {
289: /* top boundary,and i <> 0 or mx-1 */
290: tw = x[j][i-1];
291: aw = 0.5*(t0 + tw);
292: dw = PetscPowScalar(aw,beta);
293: fw = dw*(t0 - tw);
295: te = x[j][i+1];
296: ae = 0.5*(t0 + te);
297: de = PetscPowScalar(ae,beta);
298: fe = de*(te - t0);
300: ts = x[j-1][i];
301: as = 0.5*(t0 + ts);
302: ds = PetscPowScalar(as,beta);
303: fs = ds*(t0 - ts);
305: fn = zero;
307: }
309: f[j][i] = - hydhx*(fe-fw) - hxdhy*(fn-fs);
311: }
312: }
313: DMDAVecRestoreArray(da,localX,&x);
314: DMDAVecRestoreArray(da,F,&f);
315: DMRestoreLocalVector(da,&localX);
316: PetscLogFlops((22.0 + 4.0*POWFLOP)*ym*xm);
317: return(0);
318: }
319: /* -------------------- Evaluate Jacobian F(x) --------------------- */
322: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ptr)
323: {
324: AppCtx *user = (AppCtx*)ptr;
325: Mat jac = *J;
327: PetscInt i,j,mx,my,xs,ys,xm,ym;
328: PetscScalar one = 1.0,hx,hy,hxdhy,hydhx,t0,tn,ts,te,tw;
329: PetscScalar dn,ds,de,dw,an,as,ae,aw,bn,bs,be,bw,gn,gs,ge,gw;
330: PetscScalar tleft,tright,beta,bm1,coef;
331: PetscScalar v[5],**x;
332: Vec localX;
333: MatStencil col[5],row;
334: DM da;
337: SNESGetDM(snes,&da);
338: DMGetLocalVector(da,&localX);
339: *flg = SAME_NONZERO_PATTERN;
340: DMDAGetInfo(da,PETSC_NULL,&mx,&my,0,0,0,0,0,0,0,0,0,0);
341: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1);
342: hxdhy = hx/hy; hydhx = hy/hx;
343: tleft = user->tleft; tright = user->tright;
344: beta = user->beta; bm1 = user->bm1; coef = user->coef;
346: /* Get ghost points */
347: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
348: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
349: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
350: DMDAVecGetArray(da,localX,&x);
352: /* Evaluate Jacobian of function */
353: for (j=ys; j<ys+ym; j++) {
354: for (i=xs; i<xs+xm; i++) {
355: t0 = x[j][i];
357: if (i > 0 && i < mx-1 && j > 0 && j < my-1) {
359: /* general interior volume */
361: tw = x[j][i-1];
362: aw = 0.5*(t0 + tw);
363: bw = PetscPowScalar(aw,bm1);
364: /* dw = bw * aw */
365: dw = PetscPowScalar(aw,beta);
366: gw = coef*bw*(t0 - tw);
368: te = x[j][i+1];
369: ae = 0.5*(t0 + te);
370: be = PetscPowScalar(ae,bm1);
371: /* de = be * ae; */
372: de = PetscPowScalar(ae,beta);
373: ge = coef*be*(te - t0);
375: ts = x[j-1][i];
376: as = 0.5*(t0 + ts);
377: bs = PetscPowScalar(as,bm1);
378: /* ds = bs * as; */
379: ds = PetscPowScalar(as,beta);
380: gs = coef*bs*(t0 - ts);
381:
382: tn = x[j+1][i];
383: an = 0.5*(t0 + tn);
384: bn = PetscPowScalar(an,bm1);
385: /* dn = bn * an; */
386: dn = PetscPowScalar(an,beta);
387: gn = coef*bn*(tn - t0);
389: v[0] = - hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
390: v[1] = - hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
391: v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
392: v[3] = - hydhx*(de + ge); col[3].j = j; col[3].i = i+1;
393: v[4] = - hxdhy*(dn + gn); col[4].j = j+1; col[4].i = i;
394: MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
396: } else if (i == 0) {
398: /* left-hand boundary */
399: tw = tleft;
400: aw = 0.5*(t0 + tw);
401: bw = PetscPowScalar(aw,bm1);
402: /* dw = bw * aw */
403: dw = PetscPowScalar(aw,beta);
404: gw = coef*bw*(t0 - tw);
405:
406: te = x[j][i + 1];
407: ae = 0.5*(t0 + te);
408: be = PetscPowScalar(ae,bm1);
409: /* de = be * ae; */
410: de = PetscPowScalar(ae,beta);
411: ge = coef*be*(te - t0);
412:
413: /* left-hand bottom boundary */
414: if (j == 0) {
416: tn = x[j+1][i];
417: an = 0.5*(t0 + tn);
418: bn = PetscPowScalar(an,bm1);
419: /* dn = bn * an; */
420: dn = PetscPowScalar(an,beta);
421: gn = coef*bn*(tn - t0);
422:
423: v[0] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[0].j = row.j = j; col[0].i = row.i = i;
424: v[1] = - hydhx*(de + ge); col[1].j = j; col[1].i = i+1;
425: v[2] = - hxdhy*(dn + gn); col[2].j = j+1; col[2].i = i;
426: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
427:
428: /* left-hand interior boundary */
429: } else if (j < my-1) {
431: ts = x[j-1][i];
432: as = 0.5*(t0 + ts);
433: bs = PetscPowScalar(as,bm1);
434: /* ds = bs * as; */
435: ds = PetscPowScalar(as,beta);
436: gs = coef*bs*(t0 - ts);
437:
438: tn = x[j+1][i];
439: an = 0.5*(t0 + tn);
440: bn = PetscPowScalar(an,bm1);
441: /* dn = bn * an; */
442: dn = PetscPowScalar(an,beta);
443: gn = coef*bn*(tn - t0);
444:
445: v[0] = - hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
446: v[1] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
447: v[2] = - hydhx*(de + ge); col[2].j = j; col[2].i = i+1;
448: v[3] = - hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i;
449: MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
450: /* left-hand top boundary */
451: } else {
453: ts = x[j-1][i];
454: as = 0.5*(t0 + ts);
455: bs = PetscPowScalar(as,bm1);
456: /* ds = bs * as; */
457: ds = PetscPowScalar(as,beta);
458: gs = coef*bs*(t0 - ts);
459:
460: v[0] = - hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
461: v[1] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
462: v[2] = - hydhx*(de + ge); col[2].j = j; col[2].i = i+1;
463: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
464: }
466: } else if (i == mx-1) {
467:
468: /* right-hand boundary */
469: tw = x[j][i-1];
470: aw = 0.5*(t0 + tw);
471: bw = PetscPowScalar(aw,bm1);
472: /* dw = bw * aw */
473: dw = PetscPowScalar(aw,beta);
474: gw = coef*bw*(t0 - tw);
475:
476: te = tright;
477: ae = 0.5*(t0 + te);
478: be = PetscPowScalar(ae,bm1);
479: /* de = be * ae; */
480: de = PetscPowScalar(ae,beta);
481: ge = coef*be*(te - t0);
482:
483: /* right-hand bottom boundary */
484: if (j == 0) {
486: tn = x[j+1][i];
487: an = 0.5*(t0 + tn);
488: bn = PetscPowScalar(an,bm1);
489: /* dn = bn * an; */
490: dn = PetscPowScalar(an,beta);
491: gn = coef*bn*(tn - t0);
492:
493: v[0] = - hydhx*(dw - gw); col[0].j = j; col[0].i = i-1;
494: v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
495: v[2] = - hxdhy*(dn + gn); col[2].j = j+1; col[2].i = i;
496: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
497:
498: /* right-hand interior boundary */
499: } else if (j < my-1) {
501: ts = x[j-1][i];
502: as = 0.5*(t0 + ts);
503: bs = PetscPowScalar(as,bm1);
504: /* ds = bs * as; */
505: ds = PetscPowScalar(as,beta);
506: gs = coef*bs*(t0 - ts);
507:
508: tn = x[j+1][i];
509: an = 0.5*(t0 + tn);
510: bn = PetscPowScalar(an,bm1);
511: /* dn = bn * an; */
512: dn = PetscPowScalar(an,beta);
513: gn = coef*bn*(tn - t0);
514:
515: v[0] = - hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
516: v[1] = - hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
517: v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
518: v[3] = - hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i;
519: MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
520: /* right-hand top boundary */
521: } else {
523: ts = x[j-1][i];
524: as = 0.5*(t0 + ts);
525: bs = PetscPowScalar(as,bm1);
526: /* ds = bs * as; */
527: ds = PetscPowScalar(as,beta);
528: gs = coef*bs*(t0 - ts);
529:
530: v[0] = - hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
531: v[1] = - hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
532: v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
533: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
534: }
536: /* bottom boundary,and i <> 0 or mx-1 */
537: } else if (j == 0) {
539: tw = x[j][i-1];
540: aw = 0.5*(t0 + tw);
541: bw = PetscPowScalar(aw,bm1);
542: /* dw = bw * aw */
543: dw = PetscPowScalar(aw,beta);
544: gw = coef*bw*(t0 - tw);
546: te = x[j][i+1];
547: ae = 0.5*(t0 + te);
548: be = PetscPowScalar(ae,bm1);
549: /* de = be * ae; */
550: de = PetscPowScalar(ae,beta);
551: ge = coef*be*(te - t0);
553: tn = x[j+1][i];
554: an = 0.5*(t0 + tn);
555: bn = PetscPowScalar(an,bm1);
556: /* dn = bn * an; */
557: dn = PetscPowScalar(an,beta);
558: gn = coef*bn*(tn - t0);
559:
560: v[0] = - hydhx*(dw - gw); col[0].j = j; col[0].i = i-1;
561: v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
562: v[2] = - hydhx*(de + ge); col[2].j = j; col[2].i = i+1;
563: v[3] = - hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i;
564: MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
565:
566: /* top boundary,and i <> 0 or mx-1 */
567: } else if (j == my-1) {
568:
569: tw = x[j][i-1];
570: aw = 0.5*(t0 + tw);
571: bw = PetscPowScalar(aw,bm1);
572: /* dw = bw * aw */
573: dw = PetscPowScalar(aw,beta);
574: gw = coef*bw*(t0 - tw);
576: te = x[j][i+1];
577: ae = 0.5*(t0 + te);
578: be = PetscPowScalar(ae,bm1);
579: /* de = be * ae; */
580: de = PetscPowScalar(ae,beta);
581: ge = coef*be*(te - t0);
583: ts = x[j-1][i];
584: as = 0.5*(t0 + ts);
585: bs = PetscPowScalar(as,bm1);
586: /* ds = bs * as; */
587: ds = PetscPowScalar(as,beta);
588: gs = coef*bs*(t0 - ts);
590: v[0] = - hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
591: v[1] = - hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
592: v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
593: v[3] = - hydhx*(de + ge); col[3].j = j; col[3].i = i+1;
594: MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
595:
596: }
597: }
598: }
599: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
600: DMDAVecRestoreArray(da,localX,&x);
601: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
602: DMRestoreLocalVector(da,&localX);
604: PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym);
605: return(0);
606: }