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