Actual source code: ex18.c
petsc-3.13.6 2020-09-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*/
23: /*
25: This example models the partial differential equation
27: - Div(alpha* T^beta (GRAD T)) = 0.
29: where beta = 2.5 and alpha = 1.0
31: BC: T_left = 1.0, T_right = 0.1, dT/dn_top = dTdn_bottom = 0.
33: in the unit square, which is uniformly discretized in each of x and
34: y in this simple encoding. The degrees of freedom are cell centered.
36: A finite volume approximation with the usual 5-point stencil
37: is used to discretize the boundary value problem to obtain a
38: nonlinear system of equations.
40: This code was contributed by David Keyes
42: */
44: #include <petscsnes.h>
45: #include <petscdm.h>
46: #include <petscdmda.h>
48: /* User-defined Section 1.5 Writing Application Codes with PETSc context */
50: typedef struct {
51: PetscReal tleft,tright; /* Dirichlet boundary conditions */
52: PetscReal beta,bm1,coef; /* nonlinear diffusivity parameterizations */
53: } AppCtx;
55: #define POWFLOP 5 /* assume a pow() takes five flops */
57: extern PetscErrorCode FormInitialGuess(SNES,Vec,void*);
58: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
59: 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);if (ierr) return ierr;
72: /* set problem parameters */
73: user.tleft = 1.0;
74: user.tright = 0.1;
75: user.beta = 2.5;
76: PetscOptionsGetReal(NULL,NULL,"-tleft",&user.tleft,NULL);
77: PetscOptionsGetReal(NULL,NULL,"-tright",&user.tright,NULL);
78: PetscOptionsGetReal(NULL,NULL,"-beta",&user.beta,NULL);
79: user.bm1 = user.beta - 1.0;
80: 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, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,5,5,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);
91: DMSetFromOptions(da);
92: DMSetUp(da);
93: DMSetApplicationContext(da,&user);
94: SNESSetDM(snes,(DM)da);
96: /*
97: Create the nonlinear solver, and tell it the functions to use
98: */
99: SNESSetFunction(snes,NULL,FormFunction,&user);
100: SNESSetJacobian(snes,NULL,NULL,FormJacobian,&user);
101: SNESSetFromOptions(snes);
102: SNESSetComputeInitialGuess(snes,FormInitialGuess,NULL);
104: SNESSolve(snes,NULL,NULL);
105: SNESGetIterationNumber(snes,&its);
106: SNESGetLinearSolveIterations(snes,&lits);
107: litspit = ((PetscReal)lits)/((PetscReal)its);
108: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
109: PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %D\n",lits);
110: PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / SNES = %e\n",(double)litspit);
112: DMDestroy(&da);
113: SNESDestroy(&snes);
114: PetscFinalize();
115: return ierr;
116: }
117: /* -------------------- Form initial approximation ----------------- */
118: PetscErrorCode FormInitialGuess(SNES snes,Vec X,void *ctx)
119: {
120: AppCtx *user;
121: PetscInt i,j,xs,ys,xm,ym;
123: PetscReal tleft;
124: PetscScalar **x;
125: DM da;
128: SNESGetDM(snes,&da);
129: DMGetApplicationContext(da,&user);
130: tleft = user->tleft;
131: /* Get ghost points */
132: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
133: DMDAVecGetArray(da,X,&x);
135: /* Compute initial guess */
136: for (j=ys; j<ys+ym; j++) {
137: for (i=xs; i<xs+xm; i++) {
138: x[j][i] = tleft;
139: }
140: }
141: DMDAVecRestoreArray(da,X,&x);
142: return(0);
143: }
144: /* -------------------- Evaluate Function F(x) --------------------- */
145: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
146: {
147: AppCtx *user = (AppCtx*)ptr;
149: PetscInt i,j,mx,my,xs,ys,xm,ym;
150: PetscScalar zero = 0.0,one = 1.0;
151: PetscScalar hx,hy,hxdhy,hydhx;
152: 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;
153: PetscScalar tleft,tright,beta;
154: PetscScalar **x,**f;
155: Vec localX;
156: DM da;
159: SNESGetDM(snes,&da);
160: DMGetLocalVector(da,&localX);
161: DMDAGetInfo(da,NULL,&mx,&my,0,0,0,0,0,0,0,0,0,0);
162: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1);
163: hxdhy = hx/hy; hydhx = hy/hx;
164: tleft = user->tleft; tright = user->tright;
165: beta = user->beta;
167: /* Get ghost points */
168: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
169: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
170: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
171: DMDAVecGetArray(da,localX,&x);
172: DMDAVecGetArray(da,F,&f);
174: /* Evaluate function */
175: for (j=ys; j<ys+ym; j++) {
176: for (i=xs; i<xs+xm; i++) {
177: t0 = x[j][i];
179: if (i > 0 && i < mx-1 && j > 0 && j < my-1) {
181: /* general interior volume */
183: tw = x[j][i-1];
184: aw = 0.5*(t0 + tw);
185: dw = PetscPowScalar(aw,beta);
186: fw = dw*(t0 - tw);
188: te = x[j][i+1];
189: ae = 0.5*(t0 + te);
190: de = PetscPowScalar(ae,beta);
191: fe = de*(te - t0);
193: ts = x[j-1][i];
194: as = 0.5*(t0 + ts);
195: ds = PetscPowScalar(as,beta);
196: fs = ds*(t0 - ts);
198: tn = x[j+1][i];
199: an = 0.5*(t0 + tn);
200: dn = PetscPowScalar(an,beta);
201: fn = dn*(tn - t0);
203: } else if (i == 0) {
205: /* left-hand boundary */
206: tw = tleft;
207: aw = 0.5*(t0 + tw);
208: dw = PetscPowScalar(aw,beta);
209: fw = dw*(t0 - tw);
211: te = x[j][i+1];
212: ae = 0.5*(t0 + te);
213: de = PetscPowScalar(ae,beta);
214: fe = de*(te - t0);
216: if (j > 0) {
217: ts = x[j-1][i];
218: as = 0.5*(t0 + ts);
219: ds = PetscPowScalar(as,beta);
220: fs = ds*(t0 - ts);
221: } else fs = zero;
223: if (j < my-1) {
224: tn = x[j+1][i];
225: an = 0.5*(t0 + tn);
226: dn = PetscPowScalar(an,beta);
227: fn = dn*(tn - t0);
228: } else fn = zero;
230: } else if (i == mx-1) {
232: /* right-hand boundary */
233: tw = x[j][i-1];
234: aw = 0.5*(t0 + tw);
235: dw = PetscPowScalar(aw,beta);
236: fw = dw*(t0 - tw);
238: te = tright;
239: ae = 0.5*(t0 + te);
240: de = PetscPowScalar(ae,beta);
241: fe = de*(te - t0);
243: if (j > 0) {
244: ts = x[j-1][i];
245: as = 0.5*(t0 + ts);
246: ds = PetscPowScalar(as,beta);
247: fs = ds*(t0 - ts);
248: } else fs = zero;
250: if (j < my-1) {
251: tn = x[j+1][i];
252: an = 0.5*(t0 + tn);
253: dn = PetscPowScalar(an,beta);
254: fn = dn*(tn - t0);
255: } else fn = zero;
257: } else if (j == 0) {
259: /* bottom boundary,and i <> 0 or mx-1 */
260: tw = x[j][i-1];
261: aw = 0.5*(t0 + tw);
262: dw = PetscPowScalar(aw,beta);
263: fw = dw*(t0 - tw);
265: te = x[j][i+1];
266: ae = 0.5*(t0 + te);
267: de = PetscPowScalar(ae,beta);
268: fe = de*(te - t0);
270: fs = zero;
272: tn = x[j+1][i];
273: an = 0.5*(t0 + tn);
274: dn = PetscPowScalar(an,beta);
275: fn = dn*(tn - t0);
277: } else if (j == my-1) {
279: /* top boundary,and i <> 0 or mx-1 */
280: tw = x[j][i-1];
281: aw = 0.5*(t0 + tw);
282: dw = PetscPowScalar(aw,beta);
283: fw = dw*(t0 - tw);
285: te = x[j][i+1];
286: ae = 0.5*(t0 + te);
287: de = PetscPowScalar(ae,beta);
288: fe = de*(te - t0);
290: ts = x[j-1][i];
291: as = 0.5*(t0 + ts);
292: ds = PetscPowScalar(as,beta);
293: fs = ds*(t0 - ts);
295: fn = zero;
297: }
299: f[j][i] = -hydhx*(fe-fw) - hxdhy*(fn-fs);
301: }
302: }
303: DMDAVecRestoreArray(da,localX,&x);
304: DMDAVecRestoreArray(da,F,&f);
305: DMRestoreLocalVector(da,&localX);
306: PetscLogFlops((22.0 + 4.0*POWFLOP)*ym*xm);
307: return(0);
308: }
309: /* -------------------- Evaluate Jacobian F(x) --------------------- */
310: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat jac,Mat B,void *ptr)
311: {
312: AppCtx *user = (AppCtx*)ptr;
314: PetscInt i,j,mx,my,xs,ys,xm,ym;
315: PetscScalar one = 1.0,hx,hy,hxdhy,hydhx,t0,tn,ts,te,tw;
316: PetscScalar dn,ds,de,dw,an,as,ae,aw,bn,bs,be,bw,gn,gs,ge,gw;
317: PetscScalar tleft,tright,beta,bm1,coef;
318: PetscScalar v[5],**x;
319: Vec localX;
320: MatStencil col[5],row;
321: DM da;
324: SNESGetDM(snes,&da);
325: DMGetLocalVector(da,&localX);
326: DMDAGetInfo(da,NULL,&mx,&my,0,0,0,0,0,0,0,0,0,0);
327: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1);
328: hxdhy = hx/hy; hydhx = hy/hx;
329: tleft = user->tleft; tright = user->tright;
330: beta = user->beta; bm1 = user->bm1; coef = user->coef;
332: /* Get ghost points */
333: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
334: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
335: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
336: DMDAVecGetArray(da,localX,&x);
338: /* Evaluate Jacobian of function */
339: for (j=ys; j<ys+ym; j++) {
340: for (i=xs; i<xs+xm; i++) {
341: t0 = x[j][i];
343: if (i > 0 && i < mx-1 && j > 0 && j < my-1) {
345: /* general interior volume */
347: tw = x[j][i-1];
348: aw = 0.5*(t0 + tw);
349: bw = PetscPowScalar(aw,bm1);
350: /* dw = bw * aw */
351: dw = PetscPowScalar(aw,beta);
352: gw = coef*bw*(t0 - tw);
354: te = x[j][i+1];
355: ae = 0.5*(t0 + te);
356: be = PetscPowScalar(ae,bm1);
357: /* de = be * ae; */
358: de = PetscPowScalar(ae,beta);
359: ge = coef*be*(te - t0);
361: ts = x[j-1][i];
362: as = 0.5*(t0 + ts);
363: bs = PetscPowScalar(as,bm1);
364: /* ds = bs * as; */
365: ds = PetscPowScalar(as,beta);
366: gs = coef*bs*(t0 - ts);
368: tn = x[j+1][i];
369: an = 0.5*(t0 + tn);
370: bn = PetscPowScalar(an,bm1);
371: /* dn = bn * an; */
372: dn = PetscPowScalar(an,beta);
373: gn = coef*bn*(tn - t0);
375: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
376: v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
377: v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
378: v[3] = -hydhx*(de + ge); col[3].j = j; col[3].i = i+1;
379: v[4] = -hxdhy*(dn + gn); col[4].j = j+1; col[4].i = i;
380: MatSetValuesStencil(B,1,&row,5,col,v,INSERT_VALUES);
382: } else if (i == 0) {
384: /* left-hand boundary */
385: tw = tleft;
386: aw = 0.5*(t0 + tw);
387: bw = PetscPowScalar(aw,bm1);
388: /* dw = bw * aw */
389: dw = PetscPowScalar(aw,beta);
390: gw = coef*bw*(t0 - tw);
392: te = x[j][i + 1];
393: ae = 0.5*(t0 + te);
394: be = PetscPowScalar(ae,bm1);
395: /* de = be * ae; */
396: de = PetscPowScalar(ae,beta);
397: ge = coef*be*(te - t0);
399: /* left-hand bottom boundary */
400: if (j == 0) {
402: tn = x[j+1][i];
403: an = 0.5*(t0 + tn);
404: bn = PetscPowScalar(an,bm1);
405: /* dn = bn * an; */
406: dn = PetscPowScalar(an,beta);
407: gn = coef*bn*(tn - t0);
409: v[0] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[0].j = row.j = j; col[0].i = row.i = i;
410: v[1] = -hydhx*(de + ge); col[1].j = j; col[1].i = i+1;
411: v[2] = -hxdhy*(dn + gn); col[2].j = j+1; col[2].i = i;
412: MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);
414: /* left-hand interior boundary */
415: } else if (j < my-1) {
417: ts = x[j-1][i];
418: as = 0.5*(t0 + ts);
419: bs = PetscPowScalar(as,bm1);
420: /* ds = bs * as; */
421: ds = PetscPowScalar(as,beta);
422: gs = coef*bs*(t0 - ts);
424: tn = x[j+1][i];
425: an = 0.5*(t0 + tn);
426: bn = PetscPowScalar(an,bm1);
427: /* dn = bn * an; */
428: dn = PetscPowScalar(an,beta);
429: gn = coef*bn*(tn - t0);
431: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
432: v[1] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
433: v[2] = -hydhx*(de + ge); col[2].j = j; col[2].i = i+1;
434: v[3] = -hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i;
435: MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);
436: /* left-hand top boundary */
437: } else {
439: ts = x[j-1][i];
440: as = 0.5*(t0 + ts);
441: bs = PetscPowScalar(as,bm1);
442: /* ds = bs * as; */
443: ds = PetscPowScalar(as,beta);
444: gs = coef*bs*(t0 - ts);
446: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
447: v[1] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
448: v[2] = -hydhx*(de + ge); col[2].j = j; col[2].i = i+1;
449: MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);
450: }
452: } else if (i == mx-1) {
454: /* right-hand boundary */
455: tw = x[j][i-1];
456: aw = 0.5*(t0 + tw);
457: bw = PetscPowScalar(aw,bm1);
458: /* dw = bw * aw */
459: dw = PetscPowScalar(aw,beta);
460: gw = coef*bw*(t0 - tw);
462: te = tright;
463: ae = 0.5*(t0 + te);
464: be = PetscPowScalar(ae,bm1);
465: /* de = be * ae; */
466: de = PetscPowScalar(ae,beta);
467: ge = coef*be*(te - t0);
469: /* right-hand bottom boundary */
470: if (j == 0) {
472: tn = x[j+1][i];
473: an = 0.5*(t0 + tn);
474: bn = PetscPowScalar(an,bm1);
475: /* dn = bn * an; */
476: dn = PetscPowScalar(an,beta);
477: gn = coef*bn*(tn - t0);
479: v[0] = -hydhx*(dw - gw); col[0].j = j; col[0].i = i-1;
480: v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
481: v[2] = -hxdhy*(dn + gn); col[2].j = j+1; col[2].i = i;
482: MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);
484: /* right-hand interior boundary */
485: } else if (j < my-1) {
487: ts = x[j-1][i];
488: as = 0.5*(t0 + ts);
489: bs = PetscPowScalar(as,bm1);
490: /* ds = bs * as; */
491: ds = PetscPowScalar(as,beta);
492: gs = coef*bs*(t0 - ts);
494: tn = x[j+1][i];
495: an = 0.5*(t0 + tn);
496: bn = PetscPowScalar(an,bm1);
497: /* dn = bn * an; */
498: dn = PetscPowScalar(an,beta);
499: gn = coef*bn*(tn - t0);
501: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
502: v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
503: v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
504: v[3] = -hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i;
505: MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);
506: /* right-hand top boundary */
507: } else {
509: ts = x[j-1][i];
510: as = 0.5*(t0 + ts);
511: bs = PetscPowScalar(as,bm1);
512: /* ds = bs * as; */
513: ds = PetscPowScalar(as,beta);
514: gs = coef*bs*(t0 - ts);
516: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
517: v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
518: v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
519: MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES);
520: }
522: /* bottom boundary,and i <> 0 or mx-1 */
523: } else if (j == 0) {
525: tw = x[j][i-1];
526: aw = 0.5*(t0 + tw);
527: bw = PetscPowScalar(aw,bm1);
528: /* dw = bw * aw */
529: dw = PetscPowScalar(aw,beta);
530: gw = coef*bw*(t0 - tw);
532: te = x[j][i+1];
533: ae = 0.5*(t0 + te);
534: be = PetscPowScalar(ae,bm1);
535: /* de = be * ae; */
536: de = PetscPowScalar(ae,beta);
537: ge = coef*be*(te - t0);
539: tn = x[j+1][i];
540: an = 0.5*(t0 + tn);
541: bn = PetscPowScalar(an,bm1);
542: /* dn = bn * an; */
543: dn = PetscPowScalar(an,beta);
544: gn = coef*bn*(tn - t0);
546: v[0] = -hydhx*(dw - gw); col[0].j = j; col[0].i = i-1;
547: v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
548: v[2] = -hydhx*(de + ge); col[2].j = j; col[2].i = i+1;
549: v[3] = -hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i;
550: MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);
552: /* top boundary,and i <> 0 or mx-1 */
553: } else if (j == my-1) {
555: tw = x[j][i-1];
556: aw = 0.5*(t0 + tw);
557: bw = PetscPowScalar(aw,bm1);
558: /* dw = bw * aw */
559: dw = PetscPowScalar(aw,beta);
560: gw = coef*bw*(t0 - tw);
562: te = x[j][i+1];
563: ae = 0.5*(t0 + te);
564: be = PetscPowScalar(ae,bm1);
565: /* de = be * ae; */
566: de = PetscPowScalar(ae,beta);
567: ge = coef*be*(te - t0);
569: ts = x[j-1][i];
570: as = 0.5*(t0 + ts);
571: bs = PetscPowScalar(as,bm1);
572: /* ds = bs * as; */
573: ds = PetscPowScalar(as,beta);
574: gs = coef*bs*(t0 - ts);
576: v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
577: v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
578: v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
579: v[3] = -hydhx*(de + ge); col[3].j = j; col[3].i = i+1;
580: MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES);
582: }
583: }
584: }
585: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
586: DMDAVecRestoreArray(da,localX,&x);
587: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
588: DMRestoreLocalVector(da,&localX);
589: if (jac != B) {
590: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
591: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
592: }
594: PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym);
595: return(0);
596: }
600: /*TEST
602: test:
603: args: -pc_type mg -ksp_type fgmres -da_refine 2 -pc_mg_galerkin pmat -snes_view
604: requires: !single
606: TEST*/