Actual source code: ex4.c
2: /* NOTE: THIS PROGRAM HAS NOT YET BEEN SET UP IN TUTORIAL STYLE. */
4: static char help[] ="This program demonstrates use of the SNES package to solve systems of nonlinear equations on a single processor.\n\
5: Both of these examples employ\n\
6: sparse storage of the Jacobian matrices and are taken from the MINPACK-2\n\
7: test suite. By default the Bratu (SFI - solid fuel ignition) test problem\n\
8: is solved. The command line options are:\n\
9: -cavity : Solve FDC (flow in a driven cavity) problem\n\
10: -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
11: problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81)\n\
12: problem FDC: <parameter> = Reynolds number (par > 0)\n\
13: -mx <xg>, where <xg> = number of grid points in the x-direction\n\
14: -my <yg>, where <yg> = number of grid points in the y-direction\n\n";
16: /*
17: 1) Solid Fuel Ignition (SFI) problem. This problem is modeled by
18: the partial differential equation
19:
20: -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
21:
22: with boundary conditions
23:
24: u = 0 for x = 0, x = 1, y = 0, y = 1.
25:
26: A finite difference approximation with the usual 5-point stencil
27: is used to discretize the boundary value problem to obtain a nonlinear
28: system of equations.
29:
30: 2) Flow in a Driven Cavity (FDC) problem. The problem is
31: formulated as a boundary value problem, which is discretized by
32: standard finite difference approximations to obtain a system of
33: nonlinear equations.
34: */
36: #include petscsnes.h
38: typedef struct {
39: PetscReal param; /* test problem parameter */
40: PetscInt mx; /* Discretization in x-direction */
41: PetscInt my; /* Discretization in y-direction */
42: } AppCtx;
45: FormFunction1(SNES,Vec,Vec,void*),
46: FormInitialGuess1(AppCtx*,Vec);
48: FormFunction2(SNES,Vec,Vec,void*),
49: FormInitialGuess2(AppCtx*,Vec);
50:
53: int main(int argc, char **argv)
54: {
55: SNES snes; /* SNES context */
56: const SNESType method = SNESLS; /* default nonlinear solution method */
57: Vec x, r; /* solution, residual vectors */
58: Mat J; /* Jacobian matrix */
59: AppCtx user; /* user-defined application context */
60: PetscDraw draw; /* drawing context */
61: PetscInt its, N, nfails;
63: PetscTruth cavity;
64: PetscReal bratu_lambda_max = 6.81, bratu_lambda_min = 0.;
65: PetscScalar *xvalues;
67: PetscInitialize(&argc, &argv,(char *)0,help);
68: /* PetscDrawOpenX(PETSC_COMM_WORLD,0,"Solution",300,0,300,300,&draw); */
69: PetscDrawCreate(PETSC_COMM_WORLD,0,"Solution",300,0,300,300,&draw);
70: PetscDrawSetType(draw,PETSC_DRAW_X);
72: user.mx = 4;
73: user.my = 4;
74: PetscOptionsGetInt(PETSC_NULL,"-mx",&user.mx,PETSC_NULL);
75: PetscOptionsGetInt(PETSC_NULL,"-my",&user.my,PETSC_NULL);
76: PetscOptionsHasName(PETSC_NULL,"-cavity",&cavity);
77: if (cavity) user.param = 100.0;
78: else user.param = 6.0;
79: PetscOptionsGetReal(PETSC_NULL,"-par",&user.param,PETSC_NULL);
80: if (!cavity && (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min)) {
81: SETERRQ(1,"Lambda is out of range");
82: }
83: N = user.mx*user.my;
84:
85: /* Set up data structures */
86: VecCreateSeq(PETSC_COMM_SELF,N,&x);
87: VecDuplicate(x,&r);
88: MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,5,PETSC_NULL,&J);
90: /* Create nonlinear solver */
91: SNESCreate(PETSC_COMM_WORLD,&snes);
92: SNESSetType(snes,method);
94: /* Set various routines and compute initial guess for nonlinear solver */
95: if (cavity){
96: FormInitialGuess2(&user,x);
97: SNESSetFunction(snes,r,FormFunction2,(void*)&user);
98: SNESSetJacobian(snes,J,J,FormJacobian2,(void*)&user);
99: } else {
100: FormInitialGuess1(&user,x);
101: SNESSetFunction(snes,r,FormFunction1,(void*)&user);
102: SNESSetJacobian(snes,J,J,FormJacobian1,(void*)&user);
103: }
105: /* Set options and solve nonlinear system */
106: SNESSetFromOptions(snes);
107: SNESSolve(snes,PETSC_NULL,x);
108: SNESGetIterationNumber(snes,&its);
109: SNESGetNumberUnsuccessfulSteps(snes,&nfails);
111: PetscPrintf(PETSC_COMM_SELF,"number of Newton iterations = %D, ",its);
112: PetscPrintf(PETSC_COMM_SELF,"number of unsuccessful steps = %D\n\n",nfails);
113: VecGetArray(x,&xvalues);
114: PetscDrawTensorContour(draw,user.mx,user.my,0,0,xvalues);
115: VecRestoreArray(x,&xvalues);
117: /* Free data structures */
118: VecDestroy(x); VecDestroy(r);
119: MatDestroy(J); SNESDestroy(snes);
120: PetscDrawDestroy(draw);
121: PetscFinalize();
123: return 0;
124: }
125: /* ------------------------------------------------------------------ */
126: /* Bratu (Solid Fuel Ignition) Test Problem */
127: /* ------------------------------------------------------------------ */
129: /* -------------------- Form initial approximation ----------------- */
133: PetscErrorCode FormInitialGuess1(AppCtx *user,Vec X)
134: {
135: PetscInt i, j, row, mx, my;
137: PetscReal lambda, temp1, temp, hx, hy, hxdhy, hydhx,sc;
138: PetscScalar *x;
140: mx = user->mx;
141: my = user->my;
142: lambda = user->param;
144: hx = 1.0 / (PetscReal)(mx-1);
145: hy = 1.0 / (PetscReal)(my-1);
146: sc = hx*hy;
147: hxdhy = hx/hy;
148: hydhx = hy/hx;
150: VecGetArray(X,&x);
151: temp1 = lambda/(lambda + 1.0);
152: for (j=0; j<my; j++) {
153: temp = (PetscReal)(PetscMin(j,my-j-1))*hy;
154: for (i=0; i<mx; i++) {
155: row = i + j*mx;
156: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
157: x[row] = 0.0;
158: continue;
159: }
160: x[row] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));
161: }
162: }
163: VecRestoreArray(X,&x);
164: return 0;
165: }
166: /* -------------------- Evaluate Function F(x) --------------------- */
167:
170: PetscErrorCode FormFunction1(SNES snes,Vec X,Vec F,void *ptr)
171: {
172: AppCtx *user = (AppCtx*)ptr;
173: PetscInt i, j, row, mx, my;
175: PetscReal two = 2.0, one = 1.0, lambda,hx, hy, hxdhy, hydhx;
176: PetscScalar ut, ub, ul, ur, u, uxx, uyy, sc,*x,*f;
178: mx = user->mx;
179: my = user->my;
180: lambda = user->param;
182: hx = one / (PetscReal)(mx-1);
183: hy = one / (PetscReal)(my-1);
184: sc = hx*hy;
185: hxdhy = hx/hy;
186: hydhx = hy/hx;
188: VecGetArray(X,&x);
189: VecGetArray(F,&f);
190: for (j=0; j<my; j++) {
191: for (i=0; i<mx; i++) {
192: row = i + j*mx;
193: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
194: f[row] = x[row];
195: continue;
196: }
197: u = x[row];
198: ub = x[row - mx];
199: ul = x[row - 1];
200: ut = x[row + mx];
201: ur = x[row + 1];
202: uxx = (-ur + two*u - ul)*hydhx;
203: uyy = (-ut + two*u - ub)*hxdhy;
204: f[row] = uxx + uyy - sc*lambda*exp(u);
205: }
206: }
207: VecRestoreArray(X,&x);
208: VecRestoreArray(F,&f);
209: return 0;
210: }
211: /* -------------------- Evaluate Jacobian F'(x) -------------------- */
215: PetscErrorCode FormJacobian1(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
216: {
217: AppCtx *user = (AppCtx*)ptr;
218: Mat jac = *J;
219: PetscInt i, j, row, mx, my, col[5];
221: PetscScalar two = 2.0, one = 1.0, lambda, v[5],sc, *x;
222: PetscReal hx, hy, hxdhy, hydhx;
225: mx = user->mx;
226: my = user->my;
227: lambda = user->param;
229: hx = 1.0 / (PetscReal)(mx-1);
230: hy = 1.0 / (PetscReal)(my-1);
231: sc = hx*hy;
232: hxdhy = hx/hy;
233: hydhx = hy/hx;
235: VecGetArray(X,&x);
236: for (j=0; j<my; j++) {
237: for (i=0; i<mx; i++) {
238: row = i + j*mx;
239: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
240: MatSetValues(jac,1,&row,1,&row,&one,INSERT_VALUES);
241: continue;
242: }
243: v[0] = -hxdhy; col[0] = row - mx;
244: v[1] = -hydhx; col[1] = row - 1;
245: v[2] = two*(hydhx + hxdhy) - sc*lambda*exp(x[row]); col[2] = row;
246: v[3] = -hydhx; col[3] = row + 1;
247: v[4] = -hxdhy; col[4] = row + mx;
248: MatSetValues(jac,1,&row,5,col,v,INSERT_VALUES);
249: }
250: }
251: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
252: VecRestoreArray(X,&x);
253: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
254: *flag = SAME_NONZERO_PATTERN;
255: return 0;
256: }
257: /* ------------------------------------------------------------------ */
258: /* Driven Cavity Test Problem */
259: /* ------------------------------------------------------------------ */
261: /* -------------------- Form initial approximation ----------------- */
265: PetscErrorCode FormInitialGuess2(AppCtx *user,Vec X)
266: {
268: PetscInt i, j, row, mx, my;
269: PetscScalar xx,yy,*x;
270: PetscReal hx, hy;
272: mx = user->mx;
273: my = user->my;
275: /* Test for invalid input parameters */
276: if ((mx <= 0) || (my <= 0)) SETERRQ(1,0);
278: hx = 1.0 / (PetscReal)(mx-1);
279: hy = 1.0 / (PetscReal)(my-1);
281: VecGetArray(X,&x);
282: yy = 0.0;
283: for (j=0; j<my; j++) {
284: xx = 0.0;
285: for (i=0; i<mx; i++) {
286: row = i + j*mx;
287: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
288: x[row] = 0.0;
289: }
290: else {
291: x[row] = - xx*(1.0 - xx)*yy*(1.0 - yy);
292: }
293: xx = xx + hx;
294: }
295: yy = yy + hy;
296: }
297: VecRestoreArray(X,&x);
298: return 0;
299: }
300: /* -------------------- Evaluate Function F(x) --------------------- */
304: PetscErrorCode FormFunction2(SNES snes,Vec X,Vec F,void *pptr)
305: {
306: AppCtx *user = (AppCtx*)pptr;
307: PetscInt i, j, row, mx, my;
309: PetscScalar two = 2.0, zero = 0.0, pb, pbb,pbr, pl,pll,p,pr,prr;
310: PetscScalar ptl,pt,ptt,dpdy,dpdx,pblap,ptlap,rey,pbl,ptr,pllap,plap,prlap;
311: PetscScalar *x,*f, hx2, hy2, hxhy2;
312: PetscReal hx, hy;
314: mx = user->mx;
315: my = user->my;
316: hx = 1.0 / (PetscReal)(mx-1);
317: hy = 1.0 / (PetscReal)(my-1);
318: hx2 = hx*hx;
319: hy2 = hy*hy;
320: hxhy2 = hx2*hy2;
321: rey = user->param;
323: VecGetArray(X,&x);
324: VecGetArray(F,&f);
325: for (j=0; j<my; j++) {
326: for (i=0; i<mx; i++) {
327: row = i + j*mx;
328: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
329: f[row] = x[row];
330: continue;
331: }
332: if (i == 1 || j == 1) {
333: pbl = zero;
334: }
335: else {
336: pbl = x[row-mx-1];
337: }
338: if (j == 1) {
339: pb = zero;
340: pbb = x[row];
341: }
342: else if (j == 2) {
343: pb = x[row-mx];
344: pbb = zero;
345: }
346: else {
347: pb = x[row-mx];
348: pbb = x[row-2*mx];
349: }
350: if (j == 1 || i == mx-2) {
351: pbr = zero;
352: }
353: else {
354: pbr = x[row-mx+1];
355: }
356: if (i == 1) {
357: pl = zero;
358: pll = x[row];
359: }
360: else if (i == 2) {
361: pl = x[row-1];
362: pll = zero;
363: }
364: else {
365: pl = x[row-1];
366: pll = x[row-2];
367: }
368: p = x[row];
369: if (i == mx-3) {
370: pr = x[row+1];
371: prr = zero;
372: }
373: else if (i == mx-2) {
374: pr = zero;
375: prr = x[row];
376: }
377: else {
378: pr = x[row+1];
379: prr = x[row+2];
380: }
381: if (j == my-2 || i == 1) {
382: ptl = zero;
383: }
384: else {
385: ptl = x[row+mx-1];
386: }
387: if (j == my-3) {
388: pt = x[row+mx];
389: ptt = zero;
390: }
391: else if (j == my-2) {
392: pt = zero;
393: ptt = x[row] + two*hy;
394: }
395: else {
396: pt = x[row+mx];
397: ptt = x[row+2*mx];
398: }
399: if (j == my-2 || i == mx-2) {
400: ptr = zero;
401: }
402: else {
403: ptr = x[row+mx+1];
404: }
406: dpdy = (pt - pb)/(two*hy);
407: dpdx = (pr - pl)/(two*hx);
409: /* Laplacians at each point in the 5 point stencil */
410: pblap = (pbr - two*pb + pbl)/hx2 + (p - two*pb + pbb)/hy2;
411: pllap = (p - two*pl + pll)/hx2 + (ptl - two*pl + pbl)/hy2;
412: plap = (pr - two*p + pl)/hx2 + (pt - two*p + pb)/hy2;
413: prlap = (prr - two*pr + p)/hx2 + (ptr - two*pr + pbr)/hy2;
414: ptlap = (ptr - two*pt + ptl)/hx2 + (ptt - two*pt + p)/hy2;
416: f[row] = hxhy2*((prlap - two*plap + pllap)/hx2
417: + (ptlap - two*plap + pblap)/hy2
418: - rey*(dpdy*(prlap - pllap)/(two*hx)
419: - dpdx*(ptlap - pblap)/(two*hy)));
420: }
421: }
422: VecRestoreArray(X,&x);
423: VecRestoreArray(F,&f);
424: return 0;
425: }
426: /* -------------------- Evaluate Jacobian F'(x) -------------------- */
430: PetscErrorCode FormJacobian2(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *pptr)
431: {
432: AppCtx *user = (AppCtx*)pptr;
433: PetscInt i, j, row, mx, my, col;
435: PetscScalar two = 2.0, one = 1.0, zero = 0.0, pb, pbb,pbr, pl,pll,p,pr,prr;
436: PetscScalar ptl,pt,ptt,dpdy,dpdx,pblap,ptlap,rey,pbl,ptr,pllap,plap,prlap;
437: PetscScalar val,four = 4.0, three = 3.0,*x;
438: PetscReal hx, hy,hx2, hy2, hxhy2;
439: PetscTruth assembled;
441: mx = user->mx;
442: my = user->my;
443: hx = 1.0 / (PetscReal)(mx-1);
444: hy = 1.0 / (PetscReal)(my-1);
445: hx2 = hx*hx;
446: hy2 = hy*hy;
447: hxhy2 = hx2*hy2;
448: rey = user->param;
450: MatAssembled(*J,&assembled);
451: if (assembled) {
452: MatZeroEntries(*J);
453: }
454: VecGetArray(X,&x);
455: for (j=0; j<my; j++) {
456: for (i=0; i<mx; i++) {
457: row = i + j*mx;
458: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
459: MatSetValues(*J,1,&row,1,&row,&one,ADD_VALUES);
460: continue;
461: }
462: if (i == 1 || j == 1) {
463: pbl = zero;
464: }
465: else {
466: pbl = x[row-mx-1];
467: }
468: if (j == 1) {
469: pb = zero;
470: pbb = x[row];
471: }
472: else if (j == 2) {
473: pb = x[row-mx];
474: pbb = zero;
475: }
476: else {
477: pb = x[row-mx];
478: pbb = x[row-2*mx];
479: }
480: if (j == 1 || i == mx-2) {
481: pbr = zero;
482: }
483: else {
484: pbr = x[row-mx+1];
485: }
486: if (i == 1) {
487: pl = zero;
488: pll = x[row];
489: }
490: else if (i == 2) {
491: pl = x[row-1];
492: pll = zero;
493: }
494: else {
495: pl = x[row-1];
496: pll = x[row-2];
497: }
498: p = x[row];
499: if (i == mx-3) {
500: pr = x[row+1];
501: prr = zero;
502: }
503: else if (i == mx-2) {
504: pr = zero;
505: prr = x[row];
506: }
507: else {
508: pr = x[row+1];
509: prr = x[row+2];
510: }
511: if (j == my-2 || i == 1) {
512: ptl = zero;
513: }
514: else {
515: ptl = x[row+mx-1];
516: }
517: if (j == my-3) {
518: pt = x[row+mx];
519: ptt = zero;
520: }
521: else if (j == my-2) {
522: pt = zero;
523: ptt = x[row] + two*hy;
524: }
525: else {
526: pt = x[row+mx];
527: ptt = x[row+2*mx];
528: }
529: if (j == my-2 || i == mx-2) {
530: ptr = zero;
531: }
532: else {
533: ptr = x[row+mx+1];
534: }
536: dpdy = (pt - pb)/(two*hy);
537: dpdx = (pr - pl)/(two*hx);
539: /* Laplacians at each point in the 5 point stencil */
540: pblap = (pbr - two*pb + pbl)/hx2 + (p - two*pb + pbb)/hy2;
541: pllap = (p - two*pl + pll)/hx2 + (ptl - two*pl + pbl)/hy2;
542: plap = (pr - two*p + pl)/hx2 + (pt - two*p + pb)/hy2;
543: prlap = (prr - two*pr + p)/hx2 + (ptr - two*pr + pbr)/hy2;
544: ptlap = (ptr - two*pt + ptl)/hx2 + (ptt - two*pt + p)/hy2;
546: if (j > 2) {
547: val = hxhy2*(one/hy2/hy2 - rey*dpdx/hy2/(two*hy));
548: col = row - 2*mx;
549: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
550: }
551: if (i > 2) {
552: val = hxhy2*(one/hx2/hx2 + rey*dpdy/hx2/(two*hx));
553: col = row - 2;
554: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
555: }
556: if (i < mx-3) {
557: val = hxhy2*(one/hx2/hx2 - rey*dpdy/hx2/(two*hx));
558: col = row + 2;
559: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
560: }
561: if (j < my-3) {
562: val = hxhy2*(one/hy2/hy2 + rey*dpdx/hy2/(two*hy));
563: col = row + 2*mx;
564: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
565: }
566: if (i != 1 && j != 1) {
567: val = hxhy2*(two/hy2/hx2 + rey*(dpdy/hy2/(two*hx) - dpdx/hx2/(two*hy)));
568: col = row - mx - 1;
569: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
570: }
571: if (j != 1 && i != mx-2) {
572: val = hxhy2*(two/hy2/hx2 - rey*(dpdy/hy2/(two*hx) + dpdx/hx2/(two*hy)));
573: col = row - mx + 1;
574: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
575: }
576: if (j != my-2 && i != 1) {
577: val = hxhy2*(two/hy2/hx2 + rey*(dpdy/hy2/(two*hx) + dpdx/hx2/(two*hy)));
578: col = row + mx - 1;
579: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
580: }
581: if (j != my-2 && i != mx-2) {
582: val = hxhy2*(two/hy2/hx2 - rey*(dpdy/hy2/(two*hx) - dpdx/hx2/(two*hy)));
583: col = row + mx + 1;
584: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
585: }
586: if (j != 1) {
587: val = hxhy2*(-four*(one/hy2/hx2 + one/hy2/hy2)
588: + rey*((prlap - pllap)/(two*hx)/(two*hy)
589: + dpdx*(one/hx2 + one/hy2)/hy));
590: col = row - mx;
591: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
592: }
593: if (i != 1) {
594: val = hxhy2*(-four*(one/hy2/hx2 + one/hx2/hx2)
595: - rey*((ptlap - pblap)/(two*hx)/(two*hy)
596: + dpdy*(one/hx2 + one/hy2)/hx));
597: col = row - 1;
598: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
599: }
600: if (i != mx-2) {
601: val = hxhy2*(-four*(one/hy2/hx2 + one/hx2/hx2)
602: + rey*((ptlap - pblap)/(two*hx)/(two*hy)
603: + dpdy*(one/hx2 + one/hy2)/hx));
604: col = row + 1;
605: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
606: }
607: if (j != my-2) {
608: val = hxhy2*(-four*(one/hy2/hx2 + one/hy2/hy2)
609: - rey*((prlap - pllap)/(two*hx)/(two*hy)
610: + dpdx*(one/hx2 + one/hy2)/hy));
611: col = row + mx;
612: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
613: }
614: val = hxhy2*(two*(four/hx2/hy2 + three/hx2/hx2 + three/hy2/hy2));
615: MatSetValues(*J,1,&row,1,&row,&val,ADD_VALUES);
616: if (j == 1) {
617: val = hxhy2*(one/hy2/hy2 - rey*(dpdx/hy2/(two*hy)));
618: MatSetValues(*J,1,&row,1,&row,&val,ADD_VALUES);
619: }
620: if (i == 1) {
621: val = hxhy2*(one/hx2/hx2 + rey*(dpdy/hx2/(two*hx)));
622: MatSetValues(*J,1,&row,1,&row,&val,ADD_VALUES);
623: }
624: if (i == mx-2) {
625: val = hxhy2*(one/hx2/hx2 - rey*(dpdy/hx2/(two*hx)));
626: MatSetValues(*J,1,&row,1,&row,&val,ADD_VALUES);
627: }
628: if (j == my-2) {
629: val = hxhy2*(one/hy2/hy2 + rey*(dpdx/hy2/(two*hy)));
630: MatSetValues(*J,1,&row,1,&row,&val,ADD_VALUES);
631: }
632: }
633: }
634: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
635: VecRestoreArray(X,&x);
636: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
637: *flag = SAME_NONZERO_PATTERN;
638: return 0;
639: }