Actual source code: ex18.c
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*);
59: int main(int argc,char **argv)
60: {
61: SNES snes;
62: AppCtx user;
63: PetscInt its,lits;
64: PetscReal litspit;
65: DM da;
67: PetscInitialize(&argc,&argv,NULL,help);
69: /* set problem parameters */
70: user.tleft = 1.0;
71: user.tright = 0.1;
72: user.beta = 2.5;
73: PetscOptionsGetReal(NULL,NULL,"-tleft",&user.tleft,NULL);
74: PetscOptionsGetReal(NULL,NULL,"-tright",&user.tright,NULL);
75: PetscOptionsGetReal(NULL,NULL,"-beta",&user.beta,NULL);
76: user.bm1 = user.beta - 1.0;
77: user.coef = user.beta/2.0;
79: /*
80: Create the multilevel DM data structure
81: */
82: SNESCreate(PETSC_COMM_WORLD,&snes);
84: /*
85: Set the DMDA (grid structure) for the grids.
86: */
87: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,5,5,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);
88: DMSetFromOptions(da);
89: DMSetUp(da);
90: DMSetApplicationContext(da,&user);
91: SNESSetDM(snes,(DM)da);
93: /*
94: Create the nonlinear solver, and tell it the functions to use
95: */
96: SNESSetFunction(snes,NULL,FormFunction,&user);
97: SNESSetJacobian(snes,NULL,NULL,FormJacobian,&user);
98: SNESSetFromOptions(snes);
99: SNESSetComputeInitialGuess(snes,FormInitialGuess,NULL);
101: SNESSolve(snes,NULL,NULL);
102: SNESGetIterationNumber(snes,&its);
103: SNESGetLinearSolveIterations(snes,&lits);
104: litspit = ((PetscReal)lits)/((PetscReal)its);
105: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
106: PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %D\n",lits);
107: PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / SNES = %e\n",(double)litspit);
109: DMDestroy(&da);
110: SNESDestroy(&snes);
111: PetscFinalize();
112: return 0;
113: }
114: /* -------------------- Form initial approximation ----------------- */
115: PetscErrorCode FormInitialGuess(SNES snes,Vec X,void *ctx)
116: {
117: AppCtx *user;
118: PetscInt i,j,xs,ys,xm,ym;
119: PetscReal tleft;
120: PetscScalar **x;
121: DM da;
124: SNESGetDM(snes,&da);
125: DMGetApplicationContext(da,&user);
126: tleft = user->tleft;
127: /* Get ghost points */
128: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
129: DMDAVecGetArray(da,X,&x);
131: /* Compute initial guess */
132: for (j=ys; j<ys+ym; j++) {
133: for (i=xs; i<xs+xm; i++) {
134: x[j][i] = tleft;
135: }
136: }
137: DMDAVecRestoreArray(da,X,&x);
138: return 0;
139: }
140: /* -------------------- Evaluate Function F(x) --------------------- */
141: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
142: {
143: AppCtx *user = (AppCtx*)ptr;
144: PetscInt i,j,mx,my,xs,ys,xm,ym;
145: PetscScalar zero = 0.0,one = 1.0;
146: PetscScalar hx,hy,hxdhy,hydhx;
147: 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;
148: PetscScalar tleft,tright,beta;
149: PetscScalar **x,**f;
150: Vec localX;
151: DM da;
154: SNESGetDM(snes,&da);
155: DMGetLocalVector(da,&localX);
156: DMDAGetInfo(da,NULL,&mx,&my,0,0,0,0,0,0,0,0,0,0);
157: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1);
158: hxdhy = hx/hy; hydhx = hy/hx;
159: tleft = user->tleft; tright = user->tright;
160: beta = user->beta;
162: /* Get ghost points */
163: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
164: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
165: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
166: DMDAVecGetArray(da,localX,&x);
167: DMDAVecGetArray(da,F,&f);
169: /* Evaluate function */
170: for (j=ys; j<ys+ym; j++) {
171: for (i=xs; i<xs+xm; i++) {
172: t0 = x[j][i];
174: if (i > 0 && i < mx-1 && j > 0 && j < my-1) {
176: /* general interior volume */
178: tw = x[j][i-1];
179: aw = 0.5*(t0 + tw);
180: dw = PetscPowScalar(aw,beta);
181: fw = dw*(t0 - tw);
183: te = x[j][i+1];
184: ae = 0.5*(t0 + te);
185: de = PetscPowScalar(ae,beta);
186: fe = de*(te - t0);
188: ts = x[j-1][i];
189: as = 0.5*(t0 + ts);
190: ds = PetscPowScalar(as,beta);
191: fs = ds*(t0 - ts);
193: tn = x[j+1][i];
194: an = 0.5*(t0 + tn);
195: dn = PetscPowScalar(an,beta);
196: fn = dn*(tn - t0);
198: } else if (i == 0) {
200: /* left-hand boundary */
201: tw = tleft;
202: aw = 0.5*(t0 + tw);
203: dw = PetscPowScalar(aw,beta);
204: fw = dw*(t0 - tw);
206: te = x[j][i+1];
207: ae = 0.5*(t0 + te);
208: de = PetscPowScalar(ae,beta);
209: fe = de*(te - t0);
211: if (j > 0) {
212: ts = x[j-1][i];
213: as = 0.5*(t0 + ts);
214: ds = PetscPowScalar(as,beta);
215: fs = ds*(t0 - ts);
216: } else fs = zero;
218: if (j < my-1) {
219: tn = x[j+1][i];
220: an = 0.5*(t0 + tn);
221: dn = PetscPowScalar(an,beta);
222: fn = dn*(tn - t0);
223: } else fn = zero;
225: } else if (i == mx-1) {
227: /* right-hand boundary */
228: tw = x[j][i-1];
229: aw = 0.5*(t0 + tw);
230: dw = PetscPowScalar(aw,beta);
231: fw = dw*(t0 - tw);
233: te = tright;
234: ae = 0.5*(t0 + te);
235: de = PetscPowScalar(ae,beta);
236: fe = de*(te - t0);
238: if (j > 0) {
239: ts = x[j-1][i];
240: as = 0.5*(t0 + ts);
241: ds = PetscPowScalar(as,beta);
242: fs = ds*(t0 - ts);
243: } else fs = zero;
245: if (j < my-1) {
246: tn = x[j+1][i];
247: an = 0.5*(t0 + tn);
248: dn = PetscPowScalar(an,beta);
249: fn = dn*(tn - t0);
250: } else fn = zero;
252: } else if (j == 0) {
254: /* bottom boundary,and i <> 0 or mx-1 */
255: tw = x[j][i-1];
256: aw = 0.5*(t0 + tw);
257: dw = PetscPowScalar(aw,beta);
258: fw = dw*(t0 - tw);
260: te = x[j][i+1];
261: ae = 0.5*(t0 + te);
262: de = PetscPowScalar(ae,beta);
263: fe = de*(te - t0);
265: fs = zero;
267: tn = x[j+1][i];
268: an = 0.5*(t0 + tn);
269: dn = PetscPowScalar(an,beta);
270: fn = dn*(tn - t0);
272: } else if (j == my-1) {
274: /* top boundary,and i <> 0 or mx-1 */
275: tw = x[j][i-1];
276: aw = 0.5*(t0 + tw);
277: dw = PetscPowScalar(aw,beta);
278: fw = dw*(t0 - tw);
280: te = x[j][i+1];
281: ae = 0.5*(t0 + te);
282: de = PetscPowScalar(ae,beta);
283: fe = de*(te - t0);
285: ts = x[j-1][i];
286: as = 0.5*(t0 + ts);
287: ds = PetscPowScalar(as,beta);
288: fs = ds*(t0 - ts);
290: fn = zero;
292: }
294: f[j][i] = -hydhx*(fe-fw) - hxdhy*(fn-fs);
296: }
297: }
298: DMDAVecRestoreArray(da,localX,&x);
299: DMDAVecRestoreArray(da,F,&f);
300: DMRestoreLocalVector(da,&localX);
301: PetscLogFlops((22.0 + 4.0*POWFLOP)*ym*xm);
302: return 0;
303: }
304: /* -------------------- Evaluate Jacobian F(x) --------------------- */
305: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat jac,Mat B,void *ptr)
306: {
307: AppCtx *user = (AppCtx*)ptr;
308: PetscInt i,j,mx,my,xs,ys,xm,ym;
309: PetscScalar one = 1.0,hx,hy,hxdhy,hydhx,t0,tn,ts,te,tw;
310: PetscScalar dn,ds,de,dw,an,as,ae,aw,bn,bs,be,bw,gn,gs,ge,gw;
311: PetscScalar tleft,tright,beta,bm1,coef;
312: PetscScalar v[5],**x;
313: Vec localX;
314: MatStencil col[5],row;
315: DM da;
318: SNESGetDM(snes,&da);
319: DMGetLocalVector(da,&localX);
320: DMDAGetInfo(da,NULL,&mx,&my,0,0,0,0,0,0,0,0,0,0);
321: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1);
322: hxdhy = hx/hy; hydhx = hy/hx;
323: tleft = user->tleft; tright = user->tright;
324: beta = user->beta; bm1 = user->bm1; coef = user->coef;
326: /* Get ghost points */
327: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
328: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
329: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
330: DMDAVecGetArray(da,localX,&x);
332: /* Evaluate Jacobian of function */
333: for (j=ys; j<ys+ym; j++) {
334: for (i=xs; i<xs+xm; i++) {
335: t0 = x[j][i];
337: if (i > 0 && i < mx-1 && j > 0 && j < my-1) {
339: /* general interior volume */
341: tw = x[j][i-1];
342: aw = 0.5*(t0 + tw);
343: bw = PetscPowScalar(aw,bm1);
344: /* dw = bw * aw */
345: dw = PetscPowScalar(aw,beta);
346: gw = coef*bw*(t0 - tw);
348: te = x[j][i+1];
349: ae = 0.5*(t0 + te);
350: be = PetscPowScalar(ae,bm1);
351: /* de = be * ae; */
352: de = PetscPowScalar(ae,beta);
353: ge = coef*be*(te - t0);
355: ts = x[j-1][i];
356: as = 0.5*(t0 + ts);
357: bs = PetscPowScalar(as,bm1);
358: /* ds = bs * as; */
359: ds = PetscPowScalar(as,beta);
360: gs = coef*bs*(t0 - ts);
362: tn = x[j+1][i];
363: an = 0.5*(t0 + tn);
364: bn = PetscPowScalar(an,bm1);
365: /* dn = bn * an; */
366: dn = PetscPowScalar(an,beta);
367: gn = coef*bn*(tn - t0);
369: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
370: v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
371: v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
372: v[3] = -hydhx*(de + ge); col[3].j = j; col[3].i = i+1;
373: v[4] = -hxdhy*(dn + gn); col[4].j = j+1; col[4].i = i;
374: MatSetValuesStencil(B,1,&row,5,col,v,INSERT_VALUES);
376: } else if (i == 0) {
378: /* left-hand boundary */
379: tw = tleft;
380: aw = 0.5*(t0 + tw);
381: bw = PetscPowScalar(aw,bm1);
382: /* dw = bw * aw */
383: dw = PetscPowScalar(aw,beta);
384: gw = coef*bw*(t0 - tw);
386: te = x[j][i + 1];
387: ae = 0.5*(t0 + te);
388: be = PetscPowScalar(ae,bm1);
389: /* de = be * ae; */
390: de = PetscPowScalar(ae,beta);
391: ge = coef*be*(te - t0);
393: /* left-hand bottom boundary */
394: if (j == 0) {
396: tn = x[j+1][i];
397: an = 0.5*(t0 + tn);
398: bn = PetscPowScalar(an,bm1);
399: /* dn = bn * an; */
400: dn = PetscPowScalar(an,beta);
401: gn = coef*bn*(tn - t0);
403: v[0] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[0].j = row.j = j; col[0].i = row.i = i;
404: v[1] = -hydhx*(de + ge); col[1].j = j; col[1].i = i+1;
405: v[2] = -hxdhy*(dn + gn); col[2].j = j+1; col[2].i = i;
406: MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);
408: /* left-hand interior boundary */
409: } else if (j < my-1) {
411: ts = x[j-1][i];
412: as = 0.5*(t0 + ts);
413: bs = PetscPowScalar(as,bm1);
414: /* ds = bs * as; */
415: ds = PetscPowScalar(as,beta);
416: gs = coef*bs*(t0 - ts);
418: tn = x[j+1][i];
419: an = 0.5*(t0 + tn);
420: bn = PetscPowScalar(an,bm1);
421: /* dn = bn * an; */
422: dn = PetscPowScalar(an,beta);
423: gn = coef*bn*(tn - t0);
425: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
426: v[1] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
427: v[2] = -hydhx*(de + ge); col[2].j = j; col[2].i = i+1;
428: v[3] = -hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i;
429: MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);
430: /* left-hand top boundary */
431: } else {
433: ts = x[j-1][i];
434: as = 0.5*(t0 + ts);
435: bs = PetscPowScalar(as,bm1);
436: /* ds = bs * as; */
437: ds = PetscPowScalar(as,beta);
438: gs = coef*bs*(t0 - ts);
440: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
441: v[1] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
442: v[2] = -hydhx*(de + ge); col[2].j = j; col[2].i = i+1;
443: MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);
444: }
446: } else if (i == mx-1) {
448: /* right-hand boundary */
449: tw = x[j][i-1];
450: aw = 0.5*(t0 + tw);
451: bw = PetscPowScalar(aw,bm1);
452: /* dw = bw * aw */
453: dw = PetscPowScalar(aw,beta);
454: gw = coef*bw*(t0 - tw);
456: te = tright;
457: ae = 0.5*(t0 + te);
458: be = PetscPowScalar(ae,bm1);
459: /* de = be * ae; */
460: de = PetscPowScalar(ae,beta);
461: ge = coef*be*(te - t0);
463: /* right-hand bottom boundary */
464: if (j == 0) {
466: tn = x[j+1][i];
467: an = 0.5*(t0 + tn);
468: bn = PetscPowScalar(an,bm1);
469: /* dn = bn * an; */
470: dn = PetscPowScalar(an,beta);
471: gn = coef*bn*(tn - t0);
473: v[0] = -hydhx*(dw - gw); col[0].j = j; col[0].i = i-1;
474: v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
475: v[2] = -hxdhy*(dn + gn); col[2].j = j+1; col[2].i = i;
476: MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);
478: /* right-hand interior boundary */
479: } else if (j < my-1) {
481: ts = x[j-1][i];
482: as = 0.5*(t0 + ts);
483: bs = PetscPowScalar(as,bm1);
484: /* ds = bs * as; */
485: ds = PetscPowScalar(as,beta);
486: gs = coef*bs*(t0 - ts);
488: tn = x[j+1][i];
489: an = 0.5*(t0 + tn);
490: bn = PetscPowScalar(an,bm1);
491: /* dn = bn * an; */
492: dn = PetscPowScalar(an,beta);
493: gn = coef*bn*(tn - t0);
495: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
496: v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
497: v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
498: v[3] = -hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i;
499: MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);
500: /* right-hand top boundary */
501: } else {
503: ts = x[j-1][i];
504: as = 0.5*(t0 + ts);
505: bs = PetscPowScalar(as,bm1);
506: /* ds = bs * as; */
507: ds = PetscPowScalar(as,beta);
508: gs = coef*bs*(t0 - ts);
510: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
511: v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
512: v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
513: MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);
514: }
516: /* bottom boundary,and i <> 0 or mx-1 */
517: } else if (j == 0) {
519: tw = x[j][i-1];
520: aw = 0.5*(t0 + tw);
521: bw = PetscPowScalar(aw,bm1);
522: /* dw = bw * aw */
523: dw = PetscPowScalar(aw,beta);
524: gw = coef*bw*(t0 - tw);
526: te = x[j][i+1];
527: ae = 0.5*(t0 + te);
528: be = PetscPowScalar(ae,bm1);
529: /* de = be * ae; */
530: de = PetscPowScalar(ae,beta);
531: ge = coef*be*(te - t0);
533: tn = x[j+1][i];
534: an = 0.5*(t0 + tn);
535: bn = PetscPowScalar(an,bm1);
536: /* dn = bn * an; */
537: dn = PetscPowScalar(an,beta);
538: gn = coef*bn*(tn - t0);
540: v[0] = -hydhx*(dw - gw); col[0].j = j; col[0].i = i-1;
541: v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
542: v[2] = -hydhx*(de + ge); col[2].j = j; col[2].i = i+1;
543: v[3] = -hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i;
544: MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);
546: /* top boundary,and i <> 0 or mx-1 */
547: } else if (j == my-1) {
549: tw = x[j][i-1];
550: aw = 0.5*(t0 + tw);
551: bw = PetscPowScalar(aw,bm1);
552: /* dw = bw * aw */
553: dw = PetscPowScalar(aw,beta);
554: gw = coef*bw*(t0 - tw);
556: te = x[j][i+1];
557: ae = 0.5*(t0 + te);
558: be = PetscPowScalar(ae,bm1);
559: /* de = be * ae; */
560: de = PetscPowScalar(ae,beta);
561: ge = coef*be*(te - t0);
563: ts = x[j-1][i];
564: as = 0.5*(t0 + ts);
565: bs = PetscPowScalar(as,bm1);
566: /* ds = bs * as; */
567: ds = PetscPowScalar(as,beta);
568: gs = coef*bs*(t0 - ts);
570: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
571: v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
572: v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
573: v[3] = -hydhx*(de + ge); col[3].j = j; col[3].i = i+1;
574: MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);
576: }
577: }
578: }
579: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
580: DMDAVecRestoreArray(da,localX,&x);
581: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
582: DMRestoreLocalVector(da,&localX);
583: if (jac != B) {
584: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
585: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
586: }
588: PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym);
589: return 0;
590: }
592: /*TEST
594: test:
595: suffix: 1
596: args: -pc_type mg -ksp_type fgmres -da_refine 2 -pc_mg_galerkin pmat -snes_view
597: requires: !single
599: test:
600: suffix: 2
601: args: -pc_type mg -ksp_type fgmres -da_refine 2 -pc_mg_galerkin pmat -snes_view -snes_type newtontrdc
602: requires: !single
604: TEST*/