Actual source code: ex18.c
petsc-3.6.1 2015-08-06
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 <petscdm.h>
44: #include <petscdmda.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(SNES,Vec,void*);
56: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
57: extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,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,NULL,help);
72: /* set problem parameters */
73: user.tleft = 1.0;
74: user.tright = 0.1;
75: user.beta = 2.5;
76: PetscOptionsGetReal(NULL,"-tleft",&user.tleft,NULL);
77: PetscOptionsGetReal(NULL,"-tright",&user.tright,NULL);
78: PetscOptionsGetReal(NULL,"-beta",&user.beta,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, DM_BOUNDARY_NONE, DM_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,NULL,FormFunction,&user);
99: SNESSetJacobian(snes,NULL,NULL,FormJacobian,&user);
100: SNESSetFromOptions(snes);
101: SNESSetComputeInitialGuess(snes,FormInitialGuess,NULL);
103: SNESSolve(snes,NULL,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",(double)litspit);
111: DMDestroy(&da);
112: SNESDestroy(&snes);
113: PetscFinalize();
115: return 0;
116: }
117: /* -------------------- Form initial approximation ----------------- */
120: PetscErrorCode FormInitialGuess(SNES snes,Vec X,void *ctx)
121: {
122: AppCtx *user;
123: PetscInt i,j,xs,ys,xm,ym;
125: PetscReal tleft;
126: PetscScalar **x;
127: DM da;
130: SNESGetDM(snes,&da);
131: DMGetApplicationContext(da,&user);
132: tleft = user->tleft;
133: /* Get ghost points */
134: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
135: DMDAVecGetArray(da,X,&x);
137: /* Compute initial guess */
138: for (j=ys; j<ys+ym; j++) {
139: for (i=xs; i<xs+xm; i++) {
140: x[j][i] = tleft;
141: }
142: }
143: DMDAVecRestoreArray(da,X,&x);
144: return(0);
145: }
146: /* -------------------- Evaluate Function F(x) --------------------- */
149: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
150: {
151: AppCtx *user = (AppCtx*)ptr;
153: PetscInt i,j,mx,my,xs,ys,xm,ym;
154: PetscScalar zero = 0.0,one = 1.0;
155: PetscScalar hx,hy,hxdhy,hydhx;
156: 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;
157: PetscScalar tleft,tright,beta;
158: PetscScalar **x,**f;
159: Vec localX;
160: DM da;
163: SNESGetDM(snes,&da);
164: DMGetLocalVector(da,&localX);
165: DMDAGetInfo(da,NULL,&mx,&my,0,0,0,0,0,0,0,0,0,0);
166: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1);
167: hxdhy = hx/hy; hydhx = hy/hx;
168: tleft = user->tleft; tright = user->tright;
169: beta = user->beta;
171: /* Get ghost points */
172: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
173: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
174: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
175: DMDAVecGetArray(da,localX,&x);
176: DMDAVecGetArray(da,F,&f);
178: /* Evaluate function */
179: for (j=ys; j<ys+ym; j++) {
180: for (i=xs; i<xs+xm; i++) {
181: t0 = x[j][i];
183: if (i > 0 && i < mx-1 && j > 0 && j < my-1) {
185: /* general interior volume */
187: tw = x[j][i-1];
188: aw = 0.5*(t0 + tw);
189: dw = PetscPowScalar(aw,beta);
190: fw = dw*(t0 - tw);
192: te = x[j][i+1];
193: ae = 0.5*(t0 + te);
194: de = PetscPowScalar(ae,beta);
195: fe = de*(te - t0);
197: ts = x[j-1][i];
198: as = 0.5*(t0 + ts);
199: ds = PetscPowScalar(as,beta);
200: fs = ds*(t0 - ts);
202: tn = x[j+1][i];
203: an = 0.5*(t0 + tn);
204: dn = PetscPowScalar(an,beta);
205: fn = dn*(tn - t0);
207: } else if (i == 0) {
209: /* left-hand boundary */
210: tw = tleft;
211: aw = 0.5*(t0 + tw);
212: dw = PetscPowScalar(aw,beta);
213: fw = dw*(t0 - tw);
215: te = x[j][i+1];
216: ae = 0.5*(t0 + te);
217: de = PetscPowScalar(ae,beta);
218: fe = de*(te - t0);
220: if (j > 0) {
221: ts = x[j-1][i];
222: as = 0.5*(t0 + ts);
223: ds = PetscPowScalar(as,beta);
224: fs = ds*(t0 - ts);
225: } else fs = zero;
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 fn = zero;
234: } else if (i == mx-1) {
236: /* right-hand boundary */
237: tw = x[j][i-1];
238: aw = 0.5*(t0 + tw);
239: dw = PetscPowScalar(aw,beta);
240: fw = dw*(t0 - tw);
242: te = tright;
243: ae = 0.5*(t0 + te);
244: de = PetscPowScalar(ae,beta);
245: fe = de*(te - t0);
247: if (j > 0) {
248: ts = x[j-1][i];
249: as = 0.5*(t0 + ts);
250: ds = PetscPowScalar(as,beta);
251: fs = ds*(t0 - ts);
252: } else fs = zero;
254: if (j < my-1) {
255: tn = x[j+1][i];
256: an = 0.5*(t0 + tn);
257: dn = PetscPowScalar(an,beta);
258: fn = dn*(tn - t0);
259: } else fn = zero;
261: } else if (j == 0) {
263: /* bottom boundary,and i <> 0 or mx-1 */
264: tw = x[j][i-1];
265: aw = 0.5*(t0 + tw);
266: dw = PetscPowScalar(aw,beta);
267: fw = dw*(t0 - tw);
269: te = x[j][i+1];
270: ae = 0.5*(t0 + te);
271: de = PetscPowScalar(ae,beta);
272: fe = de*(te - t0);
274: fs = zero;
276: tn = x[j+1][i];
277: an = 0.5*(t0 + tn);
278: dn = PetscPowScalar(an,beta);
279: fn = dn*(tn - t0);
281: } else if (j == my-1) {
283: /* top boundary,and i <> 0 or mx-1 */
284: tw = x[j][i-1];
285: aw = 0.5*(t0 + tw);
286: dw = PetscPowScalar(aw,beta);
287: fw = dw*(t0 - tw);
289: te = x[j][i+1];
290: ae = 0.5*(t0 + te);
291: de = PetscPowScalar(ae,beta);
292: fe = de*(te - t0);
294: ts = x[j-1][i];
295: as = 0.5*(t0 + ts);
296: ds = PetscPowScalar(as,beta);
297: fs = ds*(t0 - ts);
299: fn = zero;
301: }
303: f[j][i] = -hydhx*(fe-fw) - hxdhy*(fn-fs);
305: }
306: }
307: DMDAVecRestoreArray(da,localX,&x);
308: DMDAVecRestoreArray(da,F,&f);
309: DMRestoreLocalVector(da,&localX);
310: PetscLogFlops((22.0 + 4.0*POWFLOP)*ym*xm);
311: return(0);
312: }
313: /* -------------------- Evaluate Jacobian F(x) --------------------- */
316: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat jac,Mat B,void *ptr)
317: {
318: AppCtx *user = (AppCtx*)ptr;
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: DMDAGetInfo(da,NULL,&mx,&my,0,0,0,0,0,0,0,0,0,0);
333: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1);
334: hxdhy = hx/hy; hydhx = hy/hx;
335: tleft = user->tleft; tright = user->tright;
336: beta = user->beta; bm1 = user->bm1; coef = user->coef;
338: /* Get ghost points */
339: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
340: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
341: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
342: DMDAVecGetArray(da,localX,&x);
344: /* Evaluate Jacobian of function */
345: for (j=ys; j<ys+ym; j++) {
346: for (i=xs; i<xs+xm; i++) {
347: t0 = x[j][i];
349: if (i > 0 && i < mx-1 && j > 0 && j < my-1) {
351: /* general interior volume */
353: tw = x[j][i-1];
354: aw = 0.5*(t0 + tw);
355: bw = PetscPowScalar(aw,bm1);
356: /* dw = bw * aw */
357: dw = PetscPowScalar(aw,beta);
358: gw = coef*bw*(t0 - tw);
360: te = x[j][i+1];
361: ae = 0.5*(t0 + te);
362: be = PetscPowScalar(ae,bm1);
363: /* de = be * ae; */
364: de = PetscPowScalar(ae,beta);
365: ge = coef*be*(te - t0);
367: ts = x[j-1][i];
368: as = 0.5*(t0 + ts);
369: bs = PetscPowScalar(as,bm1);
370: /* ds = bs * as; */
371: ds = PetscPowScalar(as,beta);
372: gs = coef*bs*(t0 - ts);
374: tn = x[j+1][i];
375: an = 0.5*(t0 + tn);
376: bn = PetscPowScalar(an,bm1);
377: /* dn = bn * an; */
378: dn = PetscPowScalar(an,beta);
379: gn = coef*bn*(tn - t0);
381: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
382: v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
383: v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
384: v[3] = -hydhx*(de + ge); col[3].j = j; col[3].i = i+1;
385: v[4] = -hxdhy*(dn + gn); col[4].j = j+1; col[4].i = i;
386: MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
388: } else if (i == 0) {
390: /* left-hand boundary */
391: tw = tleft;
392: aw = 0.5*(t0 + tw);
393: bw = PetscPowScalar(aw,bm1);
394: /* dw = bw * aw */
395: dw = PetscPowScalar(aw,beta);
396: gw = coef*bw*(t0 - tw);
398: te = x[j][i + 1];
399: ae = 0.5*(t0 + te);
400: be = PetscPowScalar(ae,bm1);
401: /* de = be * ae; */
402: de = PetscPowScalar(ae,beta);
403: ge = coef*be*(te - t0);
405: /* left-hand bottom boundary */
406: if (j == 0) {
408: tn = x[j+1][i];
409: an = 0.5*(t0 + tn);
410: bn = PetscPowScalar(an,bm1);
411: /* dn = bn * an; */
412: dn = PetscPowScalar(an,beta);
413: gn = coef*bn*(tn - t0);
415: v[0] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[0].j = row.j = j; col[0].i = row.i = i;
416: v[1] = -hydhx*(de + ge); col[1].j = j; col[1].i = i+1;
417: v[2] = -hxdhy*(dn + gn); col[2].j = j+1; col[2].i = i;
418: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
420: /* left-hand interior boundary */
421: } else if (j < my-1) {
423: ts = x[j-1][i];
424: as = 0.5*(t0 + ts);
425: bs = PetscPowScalar(as,bm1);
426: /* ds = bs * as; */
427: ds = PetscPowScalar(as,beta);
428: gs = coef*bs*(t0 - ts);
430: tn = x[j+1][i];
431: an = 0.5*(t0 + tn);
432: bn = PetscPowScalar(an,bm1);
433: /* dn = bn * an; */
434: dn = PetscPowScalar(an,beta);
435: gn = coef*bn*(tn - t0);
437: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
438: v[1] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
439: v[2] = -hydhx*(de + ge); col[2].j = j; col[2].i = i+1;
440: v[3] = -hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i;
441: MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
442: /* left-hand top boundary */
443: } else {
445: ts = x[j-1][i];
446: as = 0.5*(t0 + ts);
447: bs = PetscPowScalar(as,bm1);
448: /* ds = bs * as; */
449: ds = PetscPowScalar(as,beta);
450: gs = coef*bs*(t0 - ts);
452: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
453: v[1] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
454: v[2] = -hydhx*(de + ge); col[2].j = j; col[2].i = i+1;
455: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
456: }
458: } else if (i == mx-1) {
460: /* right-hand boundary */
461: tw = x[j][i-1];
462: aw = 0.5*(t0 + tw);
463: bw = PetscPowScalar(aw,bm1);
464: /* dw = bw * aw */
465: dw = PetscPowScalar(aw,beta);
466: gw = coef*bw*(t0 - tw);
468: te = tright;
469: ae = 0.5*(t0 + te);
470: be = PetscPowScalar(ae,bm1);
471: /* de = be * ae; */
472: de = PetscPowScalar(ae,beta);
473: ge = coef*be*(te - t0);
475: /* right-hand bottom boundary */
476: if (j == 0) {
478: tn = x[j+1][i];
479: an = 0.5*(t0 + tn);
480: bn = PetscPowScalar(an,bm1);
481: /* dn = bn * an; */
482: dn = PetscPowScalar(an,beta);
483: gn = coef*bn*(tn - t0);
485: v[0] = -hydhx*(dw - gw); col[0].j = j; col[0].i = i-1;
486: v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
487: v[2] = -hxdhy*(dn + gn); col[2].j = j+1; col[2].i = i;
488: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
490: /* right-hand interior boundary */
491: } else if (j < my-1) {
493: ts = x[j-1][i];
494: as = 0.5*(t0 + ts);
495: bs = PetscPowScalar(as,bm1);
496: /* ds = bs * as; */
497: ds = PetscPowScalar(as,beta);
498: gs = coef*bs*(t0 - ts);
500: tn = x[j+1][i];
501: an = 0.5*(t0 + tn);
502: bn = PetscPowScalar(an,bm1);
503: /* dn = bn * an; */
504: dn = PetscPowScalar(an,beta);
505: gn = coef*bn*(tn - t0);
507: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
508: v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
509: v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
510: v[3] = -hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i;
511: MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
512: /* right-hand top boundary */
513: } else {
515: ts = x[j-1][i];
516: as = 0.5*(t0 + ts);
517: bs = PetscPowScalar(as,bm1);
518: /* ds = bs * as; */
519: ds = PetscPowScalar(as,beta);
520: gs = coef*bs*(t0 - ts);
522: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
523: v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
524: v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
525: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
526: }
528: /* bottom boundary,and i <> 0 or mx-1 */
529: } else if (j == 0) {
531: tw = x[j][i-1];
532: aw = 0.5*(t0 + tw);
533: bw = PetscPowScalar(aw,bm1);
534: /* dw = bw * aw */
535: dw = PetscPowScalar(aw,beta);
536: gw = coef*bw*(t0 - tw);
538: te = x[j][i+1];
539: ae = 0.5*(t0 + te);
540: be = PetscPowScalar(ae,bm1);
541: /* de = be * ae; */
542: de = PetscPowScalar(ae,beta);
543: ge = coef*be*(te - t0);
545: tn = x[j+1][i];
546: an = 0.5*(t0 + tn);
547: bn = PetscPowScalar(an,bm1);
548: /* dn = bn * an; */
549: dn = PetscPowScalar(an,beta);
550: gn = coef*bn*(tn - t0);
552: v[0] = -hydhx*(dw - gw); col[0].j = j; col[0].i = i-1;
553: v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
554: v[2] = -hydhx*(de + ge); col[2].j = j; col[2].i = i+1;
555: v[3] = -hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i;
556: MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
558: /* top boundary,and i <> 0 or mx-1 */
559: } else if (j == my-1) {
561: tw = x[j][i-1];
562: aw = 0.5*(t0 + tw);
563: bw = PetscPowScalar(aw,bm1);
564: /* dw = bw * aw */
565: dw = PetscPowScalar(aw,beta);
566: gw = coef*bw*(t0 - tw);
568: te = x[j][i+1];
569: ae = 0.5*(t0 + te);
570: be = PetscPowScalar(ae,bm1);
571: /* de = be * ae; */
572: de = PetscPowScalar(ae,beta);
573: ge = coef*be*(te - t0);
575: ts = x[j-1][i];
576: as = 0.5*(t0 + ts);
577: bs = PetscPowScalar(as,bm1);
578: /* ds = bs * as; */
579: ds = PetscPowScalar(as,beta);
580: gs = coef*bs*(t0 - ts);
582: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
583: v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
584: v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
585: v[3] = -hydhx*(de + ge); col[3].j = j; col[3].i = i+1;
586: MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
588: }
589: }
590: }
591: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
592: DMDAVecRestoreArray(da,localX,&x);
593: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
594: DMRestoreLocalVector(da,&localX);
596: PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym);
597: return(0);
598: }