Actual source code: ex18.c
petsc-3.4.5 2014-06-29
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: /*
23: This example models the partial differential equation
25: - Div(alpha* T^beta (GRAD T)) = 0.
27: where beta = 2.5 and alpha = 1.0
29: BC: T_left = 1.0, T_right = 0.1, dT/dn_top = dTdn_bottom = 0.
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.
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
40: */
42: #include <petscsnes.h>
43: #include <petscdmda.h>
45: /* User-defined application context */
47: typedef struct {
48: PetscReal tleft,tright; /* Dirichlet boundary conditions */
49: PetscReal beta,bm1,coef; /* nonlinear diffusivity parameterizations */
50: } AppCtx;
52: #define POWFLOP 5 /* assume a pow() takes five flops */
54: extern PetscErrorCode FormInitialGuess(SNES,Vec,void*);
55: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
56: extern PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
60: int main(int argc,char **argv)
61: {
62: SNES snes;
63: AppCtx user;
65: PetscInt its,lits;
66: PetscReal litspit;
67: DM da;
69: PetscInitialize(&argc,&argv,NULL,help);
71: /* set problem parameters */
72: user.tleft = 1.0;
73: user.tright = 0.1;
74: user.beta = 2.5;
75: PetscOptionsGetReal(NULL,"-tleft",&user.tleft,NULL);
76: PetscOptionsGetReal(NULL,"-tright",&user.tright,NULL);
77: PetscOptionsGetReal(NULL,"-beta",&user.beta,NULL);
78: user.bm1 = user.beta - 1.0;
79: user.coef = user.beta/2.0;
82: /*
83: Create the multilevel DM data structure
84: */
85: SNESCreate(PETSC_COMM_WORLD,&snes);
87: /*
88: Set the DMDA (grid structure) for the grids.
89: */
90: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-5,-5,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);
91: DMSetApplicationContext(da,&user);
92: SNESSetDM(snes,(DM)da);
94: /*
95: Create the nonlinear solver, and tell it the functions to use
96: */
97: SNESSetFunction(snes,NULL,FormFunction,&user);
98: SNESSetJacobian(snes,NULL,NULL,FormJacobian,&user);
99: SNESSetFromOptions(snes);
100: SNESSetComputeInitialGuess(snes,FormInitialGuess,NULL);
102: SNESSolve(snes,NULL,NULL);
103: SNESGetIterationNumber(snes,&its);
104: SNESGetLinearSolveIterations(snes,&lits);
105: litspit = ((PetscReal)lits)/((PetscReal)its);
106: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
107: PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %D\n",lits);
108: PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / SNES = %e\n",litspit);
110: DMDestroy(&da);
111: SNESDestroy(&snes);
112: PetscFinalize();
114: return 0;
115: }
116: /* -------------------- Form initial approximation ----------------- */
119: PetscErrorCode FormInitialGuess(SNES snes,Vec X,void *ctx)
120: {
121: AppCtx *user;
122: PetscInt i,j,xs,ys,xm,ym;
124: PetscReal tleft;
125: PetscScalar **x;
126: DM da;
129: SNESGetDM(snes,&da);
130: DMGetApplicationContext(da,&user);
131: tleft = user->tleft;
132: /* Get ghost points */
133: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
134: DMDAVecGetArray(da,X,&x);
136: /* Compute initial guess */
137: for (j=ys; j<ys+ym; j++) {
138: for (i=xs; i<xs+xm; i++) {
139: x[j][i] = tleft;
140: }
141: }
142: DMDAVecRestoreArray(da,X,&x);
143: return(0);
144: }
145: /* -------------------- Evaluate Function F(x) --------------------- */
148: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
149: {
150: AppCtx *user = (AppCtx*)ptr;
152: PetscInt i,j,mx,my,xs,ys,xm,ym;
153: PetscScalar zero = 0.0,one = 1.0;
154: PetscScalar hx,hy,hxdhy,hydhx;
155: 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;
156: PetscScalar tleft,tright,beta;
157: PetscScalar **x,**f;
158: Vec localX;
159: DM da;
162: SNESGetDM(snes,&da);
163: DMGetLocalVector(da,&localX);
164: DMDAGetInfo(da,NULL,&mx,&my,0,0,0,0,0,0,0,0,0,0);
165: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1);
166: hxdhy = hx/hy; hydhx = hy/hx;
167: tleft = user->tleft; tright = user->tright;
168: beta = user->beta;
170: /* Get ghost points */
171: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
172: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
173: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
174: DMDAVecGetArray(da,localX,&x);
175: DMDAVecGetArray(da,F,&f);
177: /* Evaluate function */
178: for (j=ys; j<ys+ym; j++) {
179: for (i=xs; i<xs+xm; i++) {
180: t0 = x[j][i];
182: if (i > 0 && i < mx-1 && j > 0 && j < my-1) {
184: /* general interior volume */
186: tw = x[j][i-1];
187: aw = 0.5*(t0 + tw);
188: dw = PetscPowScalar(aw,beta);
189: fw = dw*(t0 - tw);
191: te = x[j][i+1];
192: ae = 0.5*(t0 + te);
193: de = PetscPowScalar(ae,beta);
194: fe = de*(te - t0);
196: ts = x[j-1][i];
197: as = 0.5*(t0 + ts);
198: ds = PetscPowScalar(as,beta);
199: fs = ds*(t0 - ts);
201: tn = x[j+1][i];
202: an = 0.5*(t0 + tn);
203: dn = PetscPowScalar(an,beta);
204: fn = dn*(tn - t0);
206: } else if (i == 0) {
208: /* left-hand boundary */
209: tw = tleft;
210: aw = 0.5*(t0 + tw);
211: dw = PetscPowScalar(aw,beta);
212: fw = dw*(t0 - tw);
214: te = x[j][i+1];
215: ae = 0.5*(t0 + te);
216: de = PetscPowScalar(ae,beta);
217: fe = de*(te - t0);
219: if (j > 0) {
220: ts = x[j-1][i];
221: as = 0.5*(t0 + ts);
222: ds = PetscPowScalar(as,beta);
223: fs = ds*(t0 - ts);
224: } else fs = zero;
226: if (j < my-1) {
227: tn = x[j+1][i];
228: an = 0.5*(t0 + tn);
229: dn = PetscPowScalar(an,beta);
230: fn = dn*(tn - t0);
231: } else fn = zero;
233: } else if (i == mx-1) {
235: /* right-hand boundary */
236: tw = x[j][i-1];
237: aw = 0.5*(t0 + tw);
238: dw = PetscPowScalar(aw,beta);
239: fw = dw*(t0 - tw);
241: te = tright;
242: ae = 0.5*(t0 + te);
243: de = PetscPowScalar(ae,beta);
244: fe = de*(te - t0);
246: if (j > 0) {
247: ts = x[j-1][i];
248: as = 0.5*(t0 + ts);
249: ds = PetscPowScalar(as,beta);
250: fs = ds*(t0 - ts);
251: } else fs = zero;
253: if (j < my-1) {
254: tn = x[j+1][i];
255: an = 0.5*(t0 + tn);
256: dn = PetscPowScalar(an,beta);
257: fn = dn*(tn - t0);
258: } else fn = zero;
260: } else if (j == 0) {
262: /* bottom boundary,and i <> 0 or mx-1 */
263: tw = x[j][i-1];
264: aw = 0.5*(t0 + tw);
265: dw = PetscPowScalar(aw,beta);
266: fw = dw*(t0 - tw);
268: te = x[j][i+1];
269: ae = 0.5*(t0 + te);
270: de = PetscPowScalar(ae,beta);
271: fe = de*(te - t0);
273: fs = zero;
275: tn = x[j+1][i];
276: an = 0.5*(t0 + tn);
277: dn = PetscPowScalar(an,beta);
278: fn = dn*(tn - t0);
280: } else if (j == my-1) {
282: /* top boundary,and i <> 0 or mx-1 */
283: tw = x[j][i-1];
284: aw = 0.5*(t0 + tw);
285: dw = PetscPowScalar(aw,beta);
286: fw = dw*(t0 - tw);
288: te = x[j][i+1];
289: ae = 0.5*(t0 + te);
290: de = PetscPowScalar(ae,beta);
291: fe = de*(te - t0);
293: ts = x[j-1][i];
294: as = 0.5*(t0 + ts);
295: ds = PetscPowScalar(as,beta);
296: fs = ds*(t0 - ts);
298: fn = zero;
300: }
302: f[j][i] = -hydhx*(fe-fw) - hxdhy*(fn-fs);
304: }
305: }
306: DMDAVecRestoreArray(da,localX,&x);
307: DMDAVecRestoreArray(da,F,&f);
308: DMRestoreLocalVector(da,&localX);
309: PetscLogFlops((22.0 + 4.0*POWFLOP)*ym*xm);
310: return(0);
311: }
312: /* -------------------- Evaluate Jacobian F(x) --------------------- */
315: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ptr)
316: {
317: AppCtx *user = (AppCtx*)ptr;
318: Mat jac = *J;
320: PetscInt i,j,mx,my,xs,ys,xm,ym;
321: PetscScalar one = 1.0,hx,hy,hxdhy,hydhx,t0,tn,ts,te,tw;
322: PetscScalar dn,ds,de,dw,an,as,ae,aw,bn,bs,be,bw,gn,gs,ge,gw;
323: PetscScalar tleft,tright,beta,bm1,coef;
324: PetscScalar v[5],**x;
325: Vec localX;
326: MatStencil col[5],row;
327: DM da;
330: SNESGetDM(snes,&da);
331: DMGetLocalVector(da,&localX);
332: *flg = SAME_NONZERO_PATTERN;
333: DMDAGetInfo(da,NULL,&mx,&my,0,0,0,0,0,0,0,0,0,0);
334: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1);
335: hxdhy = hx/hy; hydhx = hy/hx;
336: tleft = user->tleft; tright = user->tright;
337: beta = user->beta; bm1 = user->bm1; coef = user->coef;
339: /* Get ghost points */
340: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
341: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
342: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
343: DMDAVecGetArray(da,localX,&x);
345: /* Evaluate Jacobian of function */
346: for (j=ys; j<ys+ym; j++) {
347: for (i=xs; i<xs+xm; i++) {
348: t0 = x[j][i];
350: if (i > 0 && i < mx-1 && j > 0 && j < my-1) {
352: /* general interior volume */
354: tw = x[j][i-1];
355: aw = 0.5*(t0 + tw);
356: bw = PetscPowScalar(aw,bm1);
357: /* dw = bw * aw */
358: dw = PetscPowScalar(aw,beta);
359: gw = coef*bw*(t0 - tw);
361: te = x[j][i+1];
362: ae = 0.5*(t0 + te);
363: be = PetscPowScalar(ae,bm1);
364: /* de = be * ae; */
365: de = PetscPowScalar(ae,beta);
366: ge = coef*be*(te - t0);
368: ts = x[j-1][i];
369: as = 0.5*(t0 + ts);
370: bs = PetscPowScalar(as,bm1);
371: /* ds = bs * as; */
372: ds = PetscPowScalar(as,beta);
373: gs = coef*bs*(t0 - ts);
375: tn = x[j+1][i];
376: an = 0.5*(t0 + tn);
377: bn = PetscPowScalar(an,bm1);
378: /* dn = bn * an; */
379: dn = PetscPowScalar(an,beta);
380: gn = coef*bn*(tn - t0);
382: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
383: v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
384: v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
385: v[3] = -hydhx*(de + ge); col[3].j = j; col[3].i = i+1;
386: v[4] = -hxdhy*(dn + gn); col[4].j = j+1; col[4].i = i;
387: MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
389: } else if (i == 0) {
391: /* left-hand boundary */
392: tw = tleft;
393: aw = 0.5*(t0 + tw);
394: bw = PetscPowScalar(aw,bm1);
395: /* dw = bw * aw */
396: dw = PetscPowScalar(aw,beta);
397: gw = coef*bw*(t0 - tw);
399: te = x[j][i + 1];
400: ae = 0.5*(t0 + te);
401: be = PetscPowScalar(ae,bm1);
402: /* de = be * ae; */
403: de = PetscPowScalar(ae,beta);
404: ge = coef*be*(te - t0);
406: /* left-hand bottom boundary */
407: if (j == 0) {
409: tn = x[j+1][i];
410: an = 0.5*(t0 + tn);
411: bn = PetscPowScalar(an,bm1);
412: /* dn = bn * an; */
413: dn = PetscPowScalar(an,beta);
414: gn = coef*bn*(tn - t0);
416: v[0] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[0].j = row.j = j; col[0].i = row.i = i;
417: v[1] = -hydhx*(de + ge); col[1].j = j; col[1].i = i+1;
418: v[2] = -hxdhy*(dn + gn); col[2].j = j+1; col[2].i = i;
419: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
421: /* left-hand interior boundary */
422: } else if (j < my-1) {
424: ts = x[j-1][i];
425: as = 0.5*(t0 + ts);
426: bs = PetscPowScalar(as,bm1);
427: /* ds = bs * as; */
428: ds = PetscPowScalar(as,beta);
429: gs = coef*bs*(t0 - ts);
431: tn = x[j+1][i];
432: an = 0.5*(t0 + tn);
433: bn = PetscPowScalar(an,bm1);
434: /* dn = bn * an; */
435: dn = PetscPowScalar(an,beta);
436: gn = coef*bn*(tn - t0);
438: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
439: v[1] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
440: v[2] = -hydhx*(de + ge); col[2].j = j; col[2].i = i+1;
441: v[3] = -hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i;
442: MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
443: /* left-hand top boundary */
444: } else {
446: ts = x[j-1][i];
447: as = 0.5*(t0 + ts);
448: bs = PetscPowScalar(as,bm1);
449: /* ds = bs * as; */
450: ds = PetscPowScalar(as,beta);
451: gs = coef*bs*(t0 - ts);
453: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
454: v[1] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
455: v[2] = -hydhx*(de + ge); col[2].j = j; col[2].i = i+1;
456: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
457: }
459: } else if (i == mx-1) {
461: /* right-hand boundary */
462: tw = x[j][i-1];
463: aw = 0.5*(t0 + tw);
464: bw = PetscPowScalar(aw,bm1);
465: /* dw = bw * aw */
466: dw = PetscPowScalar(aw,beta);
467: gw = coef*bw*(t0 - tw);
469: te = tright;
470: ae = 0.5*(t0 + te);
471: be = PetscPowScalar(ae,bm1);
472: /* de = be * ae; */
473: de = PetscPowScalar(ae,beta);
474: ge = coef*be*(te - t0);
476: /* right-hand bottom boundary */
477: if (j == 0) {
479: tn = x[j+1][i];
480: an = 0.5*(t0 + tn);
481: bn = PetscPowScalar(an,bm1);
482: /* dn = bn * an; */
483: dn = PetscPowScalar(an,beta);
484: gn = coef*bn*(tn - t0);
486: v[0] = -hydhx*(dw - gw); col[0].j = j; col[0].i = i-1;
487: v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
488: v[2] = -hxdhy*(dn + gn); col[2].j = j+1; col[2].i = i;
489: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
491: /* right-hand interior boundary */
492: } else if (j < my-1) {
494: ts = x[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[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: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
509: v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
510: v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
511: v[3] = -hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i;
512: MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
513: /* right-hand top boundary */
514: } else {
516: ts = x[j-1][i];
517: as = 0.5*(t0 + ts);
518: bs = PetscPowScalar(as,bm1);
519: /* ds = bs * as; */
520: ds = PetscPowScalar(as,beta);
521: gs = coef*bs*(t0 - ts);
523: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
524: v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
525: v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
526: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
527: }
529: /* bottom boundary,and i <> 0 or mx-1 */
530: } else if (j == 0) {
532: tw = x[j][i-1];
533: aw = 0.5*(t0 + tw);
534: bw = PetscPowScalar(aw,bm1);
535: /* dw = bw * aw */
536: dw = PetscPowScalar(aw,beta);
537: gw = coef*bw*(t0 - tw);
539: te = x[j][i+1];
540: ae = 0.5*(t0 + te);
541: be = PetscPowScalar(ae,bm1);
542: /* de = be * ae; */
543: de = PetscPowScalar(ae,beta);
544: ge = coef*be*(te - t0);
546: tn = x[j+1][i];
547: an = 0.5*(t0 + tn);
548: bn = PetscPowScalar(an,bm1);
549: /* dn = bn * an; */
550: dn = PetscPowScalar(an,beta);
551: gn = coef*bn*(tn - t0);
553: v[0] = -hydhx*(dw - gw); col[0].j = j; col[0].i = i-1;
554: v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
555: v[2] = -hydhx*(de + ge); col[2].j = j; col[2].i = i+1;
556: v[3] = -hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i;
557: MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
559: /* top boundary,and i <> 0 or mx-1 */
560: } else if (j == my-1) {
562: tw = x[j][i-1];
563: aw = 0.5*(t0 + tw);
564: bw = PetscPowScalar(aw,bm1);
565: /* dw = bw * aw */
566: dw = PetscPowScalar(aw,beta);
567: gw = coef*bw*(t0 - tw);
569: te = x[j][i+1];
570: ae = 0.5*(t0 + te);
571: be = PetscPowScalar(ae,bm1);
572: /* de = be * ae; */
573: de = PetscPowScalar(ae,beta);
574: ge = coef*be*(te - t0);
576: ts = x[j-1][i];
577: as = 0.5*(t0 + ts);
578: bs = PetscPowScalar(as,bm1);
579: /* ds = bs * as; */
580: ds = PetscPowScalar(as,beta);
581: gs = coef*bs*(t0 - ts);
583: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
584: v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
585: v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
586: v[3] = -hydhx*(de + ge); col[3].j = j; col[3].i = i+1;
587: MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
589: }
590: }
591: }
592: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
593: DMDAVecRestoreArray(da,localX,&x);
594: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
595: DMRestoreLocalVector(da,&localX);
597: PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym);
598: return(0);
599: }