Actual source code: ex18.c
petsc-3.8.4 2018-03-24
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;
64: PetscInt its,lits;
65: PetscReal litspit;
66: DM da;
68: PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
70: /* set problem parameters */
71: user.tleft = 1.0;
72: user.tright = 0.1;
73: user.beta = 2.5;
74: PetscOptionsGetReal(NULL,NULL,"-tleft",&user.tleft,NULL);
75: PetscOptionsGetReal(NULL,NULL,"-tright",&user.tright,NULL);
76: PetscOptionsGetReal(NULL,NULL,"-beta",&user.beta,NULL);
77: user.bm1 = user.beta - 1.0;
78: user.coef = user.beta/2.0;
80: /*
81: Create the multilevel DM data structure
82: */
83: SNESCreate(PETSC_COMM_WORLD,&snes);
85: /*
86: Set the DMDA (grid structure) for the grids.
87: */
88: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,5,5,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);
89: DMSetFromOptions(da);
90: DMSetUp(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",(double)litspit);
110: DMDestroy(&da);
111: SNESDestroy(&snes);
112: PetscFinalize();
113: return ierr;
114: }
115: /* -------------------- Form initial approximation ----------------- */
116: PetscErrorCode FormInitialGuess(SNES snes,Vec X,void *ctx)
117: {
118: AppCtx *user;
119: PetscInt i,j,xs,ys,xm,ym;
121: PetscReal tleft;
122: PetscScalar **x;
123: DM da;
126: SNESGetDM(snes,&da);
127: DMGetApplicationContext(da,&user);
128: tleft = user->tleft;
129: /* Get ghost points */
130: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
131: DMDAVecGetArray(da,X,&x);
133: /* Compute initial guess */
134: for (j=ys; j<ys+ym; j++) {
135: for (i=xs; i<xs+xm; i++) {
136: x[j][i] = tleft;
137: }
138: }
139: DMDAVecRestoreArray(da,X,&x);
140: return(0);
141: }
142: /* -------------------- Evaluate Function F(x) --------------------- */
143: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
144: {
145: AppCtx *user = (AppCtx*)ptr;
147: PetscInt i,j,mx,my,xs,ys,xm,ym;
148: PetscScalar zero = 0.0,one = 1.0;
149: PetscScalar hx,hy,hxdhy,hydhx;
150: 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;
151: PetscScalar tleft,tright,beta;
152: PetscScalar **x,**f;
153: Vec localX;
154: DM da;
157: SNESGetDM(snes,&da);
158: DMGetLocalVector(da,&localX);
159: DMDAGetInfo(da,NULL,&mx,&my,0,0,0,0,0,0,0,0,0,0);
160: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1);
161: hxdhy = hx/hy; hydhx = hy/hx;
162: tleft = user->tleft; tright = user->tright;
163: beta = user->beta;
165: /* Get ghost points */
166: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
167: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
168: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
169: DMDAVecGetArray(da,localX,&x);
170: DMDAVecGetArray(da,F,&f);
172: /* Evaluate function */
173: for (j=ys; j<ys+ym; j++) {
174: for (i=xs; i<xs+xm; i++) {
175: t0 = x[j][i];
177: if (i > 0 && i < mx-1 && j > 0 && j < my-1) {
179: /* general interior volume */
181: tw = x[j][i-1];
182: aw = 0.5*(t0 + tw);
183: dw = PetscPowScalar(aw,beta);
184: fw = dw*(t0 - tw);
186: te = x[j][i+1];
187: ae = 0.5*(t0 + te);
188: de = PetscPowScalar(ae,beta);
189: fe = de*(te - t0);
191: ts = x[j-1][i];
192: as = 0.5*(t0 + ts);
193: ds = PetscPowScalar(as,beta);
194: fs = ds*(t0 - ts);
196: tn = x[j+1][i];
197: an = 0.5*(t0 + tn);
198: dn = PetscPowScalar(an,beta);
199: fn = dn*(tn - t0);
201: } else if (i == 0) {
203: /* left-hand boundary */
204: tw = tleft;
205: aw = 0.5*(t0 + tw);
206: dw = PetscPowScalar(aw,beta);
207: fw = dw*(t0 - tw);
209: te = x[j][i+1];
210: ae = 0.5*(t0 + te);
211: de = PetscPowScalar(ae,beta);
212: fe = de*(te - t0);
214: if (j > 0) {
215: ts = x[j-1][i];
216: as = 0.5*(t0 + ts);
217: ds = PetscPowScalar(as,beta);
218: fs = ds*(t0 - ts);
219: } else fs = zero;
221: if (j < my-1) {
222: tn = x[j+1][i];
223: an = 0.5*(t0 + tn);
224: dn = PetscPowScalar(an,beta);
225: fn = dn*(tn - t0);
226: } else fn = zero;
228: } else if (i == mx-1) {
230: /* right-hand boundary */
231: tw = x[j][i-1];
232: aw = 0.5*(t0 + tw);
233: dw = PetscPowScalar(aw,beta);
234: fw = dw*(t0 - tw);
236: te = tright;
237: ae = 0.5*(t0 + te);
238: de = PetscPowScalar(ae,beta);
239: fe = de*(te - t0);
241: if (j > 0) {
242: ts = x[j-1][i];
243: as = 0.5*(t0 + ts);
244: ds = PetscPowScalar(as,beta);
245: fs = ds*(t0 - ts);
246: } else fs = zero;
248: if (j < my-1) {
249: tn = x[j+1][i];
250: an = 0.5*(t0 + tn);
251: dn = PetscPowScalar(an,beta);
252: fn = dn*(tn - t0);
253: } else fn = zero;
255: } else if (j == 0) {
257: /* bottom boundary,and i <> 0 or mx-1 */
258: tw = x[j][i-1];
259: aw = 0.5*(t0 + tw);
260: dw = PetscPowScalar(aw,beta);
261: fw = dw*(t0 - tw);
263: te = x[j][i+1];
264: ae = 0.5*(t0 + te);
265: de = PetscPowScalar(ae,beta);
266: fe = de*(te - t0);
268: fs = zero;
270: tn = x[j+1][i];
271: an = 0.5*(t0 + tn);
272: dn = PetscPowScalar(an,beta);
273: fn = dn*(tn - t0);
275: } else if (j == my-1) {
277: /* top boundary,and i <> 0 or mx-1 */
278: tw = x[j][i-1];
279: aw = 0.5*(t0 + tw);
280: dw = PetscPowScalar(aw,beta);
281: fw = dw*(t0 - tw);
283: te = x[j][i+1];
284: ae = 0.5*(t0 + te);
285: de = PetscPowScalar(ae,beta);
286: fe = de*(te - t0);
288: ts = x[j-1][i];
289: as = 0.5*(t0 + ts);
290: ds = PetscPowScalar(as,beta);
291: fs = ds*(t0 - ts);
293: fn = zero;
295: }
297: f[j][i] = -hydhx*(fe-fw) - hxdhy*(fn-fs);
299: }
300: }
301: DMDAVecRestoreArray(da,localX,&x);
302: DMDAVecRestoreArray(da,F,&f);
303: DMRestoreLocalVector(da,&localX);
304: PetscLogFlops((22.0 + 4.0*POWFLOP)*ym*xm);
305: return(0);
306: }
307: /* -------------------- Evaluate Jacobian F(x) --------------------- */
308: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat jac,Mat B,void *ptr)
309: {
310: AppCtx *user = (AppCtx*)ptr;
312: PetscInt i,j,mx,my,xs,ys,xm,ym;
313: PetscScalar one = 1.0,hx,hy,hxdhy,hydhx,t0,tn,ts,te,tw;
314: PetscScalar dn,ds,de,dw,an,as,ae,aw,bn,bs,be,bw,gn,gs,ge,gw;
315: PetscScalar tleft,tright,beta,bm1,coef;
316: PetscScalar v[5],**x;
317: Vec localX;
318: MatStencil col[5],row;
319: DM da;
322: SNESGetDM(snes,&da);
323: DMGetLocalVector(da,&localX);
324: DMDAGetInfo(da,NULL,&mx,&my,0,0,0,0,0,0,0,0,0,0);
325: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1);
326: hxdhy = hx/hy; hydhx = hy/hx;
327: tleft = user->tleft; tright = user->tright;
328: beta = user->beta; bm1 = user->bm1; coef = user->coef;
330: /* Get ghost points */
331: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
332: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
333: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
334: DMDAVecGetArray(da,localX,&x);
336: /* Evaluate Jacobian of function */
337: for (j=ys; j<ys+ym; j++) {
338: for (i=xs; i<xs+xm; i++) {
339: t0 = x[j][i];
341: if (i > 0 && i < mx-1 && j > 0 && j < my-1) {
343: /* general interior volume */
345: tw = x[j][i-1];
346: aw = 0.5*(t0 + tw);
347: bw = PetscPowScalar(aw,bm1);
348: /* dw = bw * aw */
349: dw = PetscPowScalar(aw,beta);
350: gw = coef*bw*(t0 - tw);
352: te = x[j][i+1];
353: ae = 0.5*(t0 + te);
354: be = PetscPowScalar(ae,bm1);
355: /* de = be * ae; */
356: de = PetscPowScalar(ae,beta);
357: ge = coef*be*(te - t0);
359: ts = x[j-1][i];
360: as = 0.5*(t0 + ts);
361: bs = PetscPowScalar(as,bm1);
362: /* ds = bs * as; */
363: ds = PetscPowScalar(as,beta);
364: gs = coef*bs*(t0 - ts);
366: tn = x[j+1][i];
367: an = 0.5*(t0 + tn);
368: bn = PetscPowScalar(an,bm1);
369: /* dn = bn * an; */
370: dn = PetscPowScalar(an,beta);
371: gn = coef*bn*(tn - t0);
373: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
374: v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
375: v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
376: v[3] = -hydhx*(de + ge); col[3].j = j; col[3].i = i+1;
377: v[4] = -hxdhy*(dn + gn); col[4].j = j+1; col[4].i = i;
378: MatSetValuesStencil(B,1,&row,5,col,v,INSERT_VALUES);
380: } else if (i == 0) {
382: /* left-hand boundary */
383: tw = tleft;
384: aw = 0.5*(t0 + tw);
385: bw = PetscPowScalar(aw,bm1);
386: /* dw = bw * aw */
387: dw = PetscPowScalar(aw,beta);
388: gw = coef*bw*(t0 - tw);
390: te = x[j][i + 1];
391: ae = 0.5*(t0 + te);
392: be = PetscPowScalar(ae,bm1);
393: /* de = be * ae; */
394: de = PetscPowScalar(ae,beta);
395: ge = coef*be*(te - t0);
397: /* left-hand bottom boundary */
398: if (j == 0) {
400: tn = x[j+1][i];
401: an = 0.5*(t0 + tn);
402: bn = PetscPowScalar(an,bm1);
403: /* dn = bn * an; */
404: dn = PetscPowScalar(an,beta);
405: gn = coef*bn*(tn - t0);
407: v[0] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[0].j = row.j = j; col[0].i = row.i = i;
408: v[1] = -hydhx*(de + ge); col[1].j = j; col[1].i = i+1;
409: v[2] = -hxdhy*(dn + gn); col[2].j = j+1; col[2].i = i;
410: MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);
412: /* left-hand interior boundary */
413: } else if (j < my-1) {
415: ts = x[j-1][i];
416: as = 0.5*(t0 + ts);
417: bs = PetscPowScalar(as,bm1);
418: /* ds = bs * as; */
419: ds = PetscPowScalar(as,beta);
420: gs = coef*bs*(t0 - ts);
422: tn = x[j+1][i];
423: an = 0.5*(t0 + tn);
424: bn = PetscPowScalar(an,bm1);
425: /* dn = bn * an; */
426: dn = PetscPowScalar(an,beta);
427: gn = coef*bn*(tn - t0);
429: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
430: v[1] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
431: v[2] = -hydhx*(de + ge); col[2].j = j; col[2].i = i+1;
432: v[3] = -hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i;
433: MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);
434: /* left-hand top boundary */
435: } else {
437: ts = x[j-1][i];
438: as = 0.5*(t0 + ts);
439: bs = PetscPowScalar(as,bm1);
440: /* ds = bs * as; */
441: ds = PetscPowScalar(as,beta);
442: gs = coef*bs*(t0 - ts);
444: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
445: v[1] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
446: v[2] = -hydhx*(de + ge); col[2].j = j; col[2].i = i+1;
447: MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);
448: }
450: } else if (i == mx-1) {
452: /* right-hand boundary */
453: tw = x[j][i-1];
454: aw = 0.5*(t0 + tw);
455: bw = PetscPowScalar(aw,bm1);
456: /* dw = bw * aw */
457: dw = PetscPowScalar(aw,beta);
458: gw = coef*bw*(t0 - tw);
460: te = tright;
461: ae = 0.5*(t0 + te);
462: be = PetscPowScalar(ae,bm1);
463: /* de = be * ae; */
464: de = PetscPowScalar(ae,beta);
465: ge = coef*be*(te - t0);
467: /* right-hand bottom boundary */
468: if (j == 0) {
470: tn = x[j+1][i];
471: an = 0.5*(t0 + tn);
472: bn = PetscPowScalar(an,bm1);
473: /* dn = bn * an; */
474: dn = PetscPowScalar(an,beta);
475: gn = coef*bn*(tn - t0);
477: v[0] = -hydhx*(dw - gw); col[0].j = j; col[0].i = i-1;
478: v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
479: v[2] = -hxdhy*(dn + gn); col[2].j = j+1; col[2].i = i;
480: MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);
482: /* right-hand interior boundary */
483: } else if (j < my-1) {
485: ts = x[j-1][i];
486: as = 0.5*(t0 + ts);
487: bs = PetscPowScalar(as,bm1);
488: /* ds = bs * as; */
489: ds = PetscPowScalar(as,beta);
490: gs = coef*bs*(t0 - ts);
492: tn = x[j+1][i];
493: an = 0.5*(t0 + tn);
494: bn = PetscPowScalar(an,bm1);
495: /* dn = bn * an; */
496: dn = PetscPowScalar(an,beta);
497: gn = coef*bn*(tn - t0);
499: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
500: v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
501: v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
502: v[3] = -hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i;
503: MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);
504: /* right-hand top boundary */
505: } else {
507: ts = x[j-1][i];
508: as = 0.5*(t0 + ts);
509: bs = PetscPowScalar(as,bm1);
510: /* ds = bs * as; */
511: ds = PetscPowScalar(as,beta);
512: gs = coef*bs*(t0 - ts);
514: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
515: v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
516: v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
517: MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);
518: }
520: /* bottom boundary,and i <> 0 or mx-1 */
521: } else if (j == 0) {
523: tw = x[j][i-1];
524: aw = 0.5*(t0 + tw);
525: bw = PetscPowScalar(aw,bm1);
526: /* dw = bw * aw */
527: dw = PetscPowScalar(aw,beta);
528: gw = coef*bw*(t0 - tw);
530: te = x[j][i+1];
531: ae = 0.5*(t0 + te);
532: be = PetscPowScalar(ae,bm1);
533: /* de = be * ae; */
534: de = PetscPowScalar(ae,beta);
535: ge = coef*be*(te - t0);
537: tn = x[j+1][i];
538: an = 0.5*(t0 + tn);
539: bn = PetscPowScalar(an,bm1);
540: /* dn = bn * an; */
541: dn = PetscPowScalar(an,beta);
542: gn = coef*bn*(tn - t0);
544: v[0] = -hydhx*(dw - gw); col[0].j = j; col[0].i = i-1;
545: v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
546: v[2] = -hydhx*(de + ge); col[2].j = j; col[2].i = i+1;
547: v[3] = -hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i;
548: MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);
550: /* top boundary,and i <> 0 or mx-1 */
551: } else if (j == my-1) {
553: tw = x[j][i-1];
554: aw = 0.5*(t0 + tw);
555: bw = PetscPowScalar(aw,bm1);
556: /* dw = bw * aw */
557: dw = PetscPowScalar(aw,beta);
558: gw = coef*bw*(t0 - tw);
560: te = x[j][i+1];
561: ae = 0.5*(t0 + te);
562: be = PetscPowScalar(ae,bm1);
563: /* de = be * ae; */
564: de = PetscPowScalar(ae,beta);
565: ge = coef*be*(te - t0);
567: ts = x[j-1][i];
568: as = 0.5*(t0 + ts);
569: bs = PetscPowScalar(as,bm1);
570: /* ds = bs * as; */
571: ds = PetscPowScalar(as,beta);
572: gs = coef*bs*(t0 - ts);
574: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
575: v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
576: v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
577: v[3] = -hydhx*(de + ge); col[3].j = j; col[3].i = i+1;
578: MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);
580: }
581: }
582: }
583: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
584: DMDAVecRestoreArray(da,localX,&x);
585: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
586: DMRestoreLocalVector(da,&localX);
587: if (jac != B) {
588: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
589: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
590: }
592: PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym);
593: return(0);
594: }