Actual source code: ex4.c
1: /*$Id: ex4.c,v 1.64 2001/08/07 21:31:12 bsmith Exp $*/
3: /* NOTE: THIS PROGRAM HAS NOT YET BEEN SET UP IN TUTORIAL STYLE. */
5: static char help[] ="This program demonstrates use of the SNES package to solve systems of nonlinear equations on a single processor.\n\
6: Both of these examples employ\n\
7: sparse storage of the Jacobian matrices and are taken from the MINPACK-2\n\
8: test suite. By default the Bratu (SFI - solid fuel ignition) test problem\n\
9: is solved. The command line options are:\n\
10: -cavity : Solve FDC (flow in a driven cavity) problem\n\
11: -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
12: problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81)\n\
13: problem FDC: <parameter> = Reynolds number (par > 0)\n\
14: -mx <xg>, where <xg> = number of grid points in the x-direction\n\
15: -my <yg>, where <yg> = number of grid points in the y-direction\n\n";
17: /*
18: 1) Solid Fuel Ignition (SFI) problem. This problem is modeled by
19: the partial differential equation
20:
21: -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
22:
23: with boundary conditions
24:
25: u = 0 for x = 0, x = 1, y = 0, y = 1.
26:
27: A finite difference approximation with the usual 5-point stencil
28: is used to discretize the boundary value problem to obtain a nonlinear
29: system of equations.
30:
31: 2) Flow in a Driven Cavity (FDC) problem. The problem is
32: formulated as a boundary value problem, which is discretized by
33: standard finite difference approximations to obtain a system of
34: nonlinear equations.
35: */
37: #include petscsnes.h
39: typedef struct {
40: PetscReal param; /* test problem parameter */
41: int mx; /* Discretization in x-direction */
42: int my; /* Discretization in y-direction */
43: } AppCtx;
45: extern int FormJacobian1(SNES,Vec,Mat*,Mat*,MatStructure*,void*),
46: FormFunction1(SNES,Vec,Vec,void*),
47: FormInitialGuess1(AppCtx*,Vec);
48: extern int FormJacobian2(SNES,Vec,Mat*,Mat*,MatStructure*,void*),
49: FormFunction2(SNES,Vec,Vec,void*),
50: FormInitialGuess2(AppCtx*,Vec);
51:
54: int main(int argc, char **argv)
55: {
56: SNES snes; /* SNES context */
57: SNESType method = SNESLS; /* default nonlinear solution method */
58: Vec x, r; /* solution, residual vectors */
59: Mat J; /* Jacobian matrix */
60: AppCtx user; /* user-defined application context */
61: PetscDraw draw; /* drawing context */
62: int ierr, its, N, nfails;
63: PetscTruth flg,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,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: int FormInitialGuess1(AppCtx *user,Vec X)
134: {
135: int i, j, row, mx, my, ierr;
136: PetscReal lambda, temp1, temp, hx, hy, hxdhy, hydhx,sc;
137: PetscScalar *x;
139: mx = user->mx;
140: my = user->my;
141: lambda = user->param;
143: hx = 1.0 / (PetscReal)(mx-1);
144: hy = 1.0 / (PetscReal)(my-1);
145: sc = hx*hy;
146: hxdhy = hx/hy;
147: hydhx = hy/hx;
149: VecGetArray(X,&x);
150: temp1 = lambda/(lambda + 1.0);
151: for (j=0; j<my; j++) {
152: temp = (PetscReal)(PetscMin(j,my-j-1))*hy;
153: for (i=0; i<mx; i++) {
154: row = i + j*mx;
155: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
156: x[row] = 0.0;
157: continue;
158: }
159: x[row] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));
160: }
161: }
162: VecRestoreArray(X,&x);
163: return 0;
164: }
165: /* -------------------- Evaluate Function F(x) --------------------- */
166:
169: int FormFunction1(SNES snes,Vec X,Vec F,void *ptr)
170: {
171: AppCtx *user = (AppCtx*)ptr;
172: int ierr, i, j, row, mx, my;
173: PetscReal two = 2.0, one = 1.0, lambda,hx, hy, hxdhy, hydhx;
174: PetscScalar ut, ub, ul, ur, u, uxx, uyy, sc,*x,*f;
176: mx = user->mx;
177: my = user->my;
178: lambda = user->param;
180: hx = one / (PetscReal)(mx-1);
181: hy = one / (PetscReal)(my-1);
182: sc = hx*hy;
183: hxdhy = hx/hy;
184: hydhx = hy/hx;
186: VecGetArray(X,&x);
187: VecGetArray(F,&f);
188: for (j=0; j<my; j++) {
189: for (i=0; i<mx; i++) {
190: row = i + j*mx;
191: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
192: f[row] = x[row];
193: continue;
194: }
195: u = x[row];
196: ub = x[row - mx];
197: ul = x[row - 1];
198: ut = x[row + mx];
199: ur = x[row + 1];
200: uxx = (-ur + two*u - ul)*hydhx;
201: uyy = (-ut + two*u - ub)*hxdhy;
202: f[row] = uxx + uyy - sc*lambda*exp(u);
203: }
204: }
205: VecRestoreArray(X,&x);
206: VecRestoreArray(F,&f);
207: return 0;
208: }
209: /* -------------------- Evaluate Jacobian F'(x) -------------------- */
213: int FormJacobian1(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
214: {
215: AppCtx *user = (AppCtx*)ptr;
216: Mat jac = *J;
217: int i, j, row, mx, my, col[5], ierr;
218: PetscScalar two = 2.0, one = 1.0, lambda, v[5],sc, *x;
219: PetscReal hx, hy, hxdhy, hydhx;
222: mx = user->mx;
223: my = user->my;
224: lambda = user->param;
226: hx = 1.0 / (PetscReal)(mx-1);
227: hy = 1.0 / (PetscReal)(my-1);
228: sc = hx*hy;
229: hxdhy = hx/hy;
230: hydhx = hy/hx;
232: VecGetArray(X,&x);
233: for (j=0; j<my; j++) {
234: for (i=0; i<mx; i++) {
235: row = i + j*mx;
236: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
237: MatSetValues(jac,1,&row,1,&row,&one,INSERT_VALUES);
238: continue;
239: }
240: v[0] = -hxdhy; col[0] = row - mx;
241: v[1] = -hydhx; col[1] = row - 1;
242: v[2] = two*(hydhx + hxdhy) - sc*lambda*exp(x[row]); col[2] = row;
243: v[3] = -hydhx; col[3] = row + 1;
244: v[4] = -hxdhy; col[4] = row + mx;
245: MatSetValues(jac,1,&row,5,col,v,INSERT_VALUES);
246: }
247: }
248: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
249: VecRestoreArray(X,&x);
250: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
251: *flag = SAME_NONZERO_PATTERN;
252: return 0;
253: }
254: /* ------------------------------------------------------------------ */
255: /* Driven Cavity Test Problem */
256: /* ------------------------------------------------------------------ */
258: /* -------------------- Form initial approximation ----------------- */
262: int FormInitialGuess2(AppCtx *user,Vec X)
263: {
264: int ierr, i, j, row, mx, my;
265: PetscScalar xx,yy,*x;
266: PetscReal hx, hy;
268: mx = user->mx;
269: my = user->my;
271: /* Test for invalid input parameters */
272: if ((mx <= 0) || (my <= 0)) SETERRQ(1,0);
274: hx = 1.0 / (PetscReal)(mx-1);
275: hy = 1.0 / (PetscReal)(my-1);
277: VecGetArray(X,&x);
278: yy = 0.0;
279: for (j=0; j<my; j++) {
280: xx = 0.0;
281: for (i=0; i<mx; i++) {
282: row = i + j*mx;
283: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
284: x[row] = 0.0;
285: }
286: else {
287: x[row] = - xx*(1.0 - xx)*yy*(1.0 - yy);
288: }
289: xx = xx + hx;
290: }
291: yy = yy + hy;
292: }
293: VecRestoreArray(X,&x);
294: return 0;
295: }
296: /* -------------------- Evaluate Function F(x) --------------------- */
300: int FormFunction2(SNES snes,Vec X,Vec F,void *pptr)
301: {
302: AppCtx *user = (AppCtx*)pptr;
303: int i, j, row, mx, my, ierr;
304: PetscScalar two = 2.0, zero = 0.0, pb, pbb,pbr, pl,pll,p,pr,prr;
305: PetscScalar ptl,pt,ptt,dpdy,dpdx,pblap,ptlap,rey,pbl,ptr,pllap,plap,prlap;
306: PetscScalar *x,*f, hx2, hy2, hxhy2;
307: PetscReal hx, hy;
309: mx = user->mx;
310: my = user->my;
311: hx = 1.0 / (PetscReal)(mx-1);
312: hy = 1.0 / (PetscReal)(my-1);
313: hx2 = hx*hx;
314: hy2 = hy*hy;
315: hxhy2 = hx2*hy2;
316: rey = user->param;
318: VecGetArray(X,&x);
319: VecGetArray(F,&f);
320: for (j=0; j<my; j++) {
321: for (i=0; i<mx; i++) {
322: row = i + j*mx;
323: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
324: f[row] = x[row];
325: continue;
326: }
327: if (i == 1 || j == 1) {
328: pbl = zero;
329: }
330: else {
331: pbl = x[row-mx-1];
332: }
333: if (j == 1) {
334: pb = zero;
335: pbb = x[row];
336: }
337: else if (j == 2) {
338: pb = x[row-mx];
339: pbb = zero;
340: }
341: else {
342: pb = x[row-mx];
343: pbb = x[row-2*mx];
344: }
345: if (j == 1 || i == mx-2) {
346: pbr = zero;
347: }
348: else {
349: pbr = x[row-mx+1];
350: }
351: if (i == 1) {
352: pl = zero;
353: pll = x[row];
354: }
355: else if (i == 2) {
356: pl = x[row-1];
357: pll = zero;
358: }
359: else {
360: pl = x[row-1];
361: pll = x[row-2];
362: }
363: p = x[row];
364: if (i == mx-3) {
365: pr = x[row+1];
366: prr = zero;
367: }
368: else if (i == mx-2) {
369: pr = zero;
370: prr = x[row];
371: }
372: else {
373: pr = x[row+1];
374: prr = x[row+2];
375: }
376: if (j == my-2 || i == 1) {
377: ptl = zero;
378: }
379: else {
380: ptl = x[row+mx-1];
381: }
382: if (j == my-3) {
383: pt = x[row+mx];
384: ptt = zero;
385: }
386: else if (j == my-2) {
387: pt = zero;
388: ptt = x[row] + two*hy;
389: }
390: else {
391: pt = x[row+mx];
392: ptt = x[row+2*mx];
393: }
394: if (j == my-2 || i == mx-2) {
395: ptr = zero;
396: }
397: else {
398: ptr = x[row+mx+1];
399: }
401: dpdy = (pt - pb)/(two*hy);
402: dpdx = (pr - pl)/(two*hx);
404: /* Laplacians at each point in the 5 point stencil */
405: pblap = (pbr - two*pb + pbl)/hx2 + (p - two*pb + pbb)/hy2;
406: pllap = (p - two*pl + pll)/hx2 + (ptl - two*pl + pbl)/hy2;
407: plap = (pr - two*p + pl)/hx2 + (pt - two*p + pb)/hy2;
408: prlap = (prr - two*pr + p)/hx2 + (ptr - two*pr + pbr)/hy2;
409: ptlap = (ptr - two*pt + ptl)/hx2 + (ptt - two*pt + p)/hy2;
411: f[row] = hxhy2*((prlap - two*plap + pllap)/hx2
412: + (ptlap - two*plap + pblap)/hy2
413: - rey*(dpdy*(prlap - pllap)/(two*hx)
414: - dpdx*(ptlap - pblap)/(two*hy)));
415: }
416: }
417: VecRestoreArray(X,&x);
418: VecRestoreArray(F,&f);
419: return 0;
420: }
421: /* -------------------- Evaluate Jacobian F'(x) -------------------- */
425: int FormJacobian2(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *pptr)
426: {
427: AppCtx *user = (AppCtx*)pptr;
428: int i, j, row, mx, my, col, ierr;
429: PetscScalar two = 2.0, one = 1.0, zero = 0.0, pb, pbb,pbr, pl,pll,p,pr,prr;
430: PetscScalar ptl,pt,ptt,dpdy,dpdx,pblap,ptlap,rey,pbl,ptr,pllap,plap,prlap;
431: PetscScalar val,four = 4.0, three = 3.0,*x;
432: PetscReal hx, hy,hx2, hy2, hxhy2;
434: mx = user->mx;
435: my = user->my;
436: hx = 1.0 / (PetscReal)(mx-1);
437: hy = 1.0 / (PetscReal)(my-1);
438: hx2 = hx*hx;
439: hy2 = hy*hy;
440: hxhy2 = hx2*hy2;
441: rey = user->param;
443: MatZeroEntries(*J);
444: VecGetArray(X,&x);
445: for (j=0; j<my; j++) {
446: for (i=0; i<mx; i++) {
447: row = i + j*mx;
448: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
449: MatSetValues(*J,1,&row,1,&row,&one,ADD_VALUES);
450: continue;
451: }
452: if (i == 1 || j == 1) {
453: pbl = zero;
454: }
455: else {
456: pbl = x[row-mx-1];
457: }
458: if (j == 1) {
459: pb = zero;
460: pbb = x[row];
461: }
462: else if (j == 2) {
463: pb = x[row-mx];
464: pbb = zero;
465: }
466: else {
467: pb = x[row-mx];
468: pbb = x[row-2*mx];
469: }
470: if (j == 1 || i == mx-2) {
471: pbr = zero;
472: }
473: else {
474: pbr = x[row-mx+1];
475: }
476: if (i == 1) {
477: pl = zero;
478: pll = x[row];
479: }
480: else if (i == 2) {
481: pl = x[row-1];
482: pll = zero;
483: }
484: else {
485: pl = x[row-1];
486: pll = x[row-2];
487: }
488: p = x[row];
489: if (i == mx-3) {
490: pr = x[row+1];
491: prr = zero;
492: }
493: else if (i == mx-2) {
494: pr = zero;
495: prr = x[row];
496: }
497: else {
498: pr = x[row+1];
499: prr = x[row+2];
500: }
501: if (j == my-2 || i == 1) {
502: ptl = zero;
503: }
504: else {
505: ptl = x[row+mx-1];
506: }
507: if (j == my-3) {
508: pt = x[row+mx];
509: ptt = zero;
510: }
511: else if (j == my-2) {
512: pt = zero;
513: ptt = x[row] + two*hy;
514: }
515: else {
516: pt = x[row+mx];
517: ptt = x[row+2*mx];
518: }
519: if (j == my-2 || i == mx-2) {
520: ptr = zero;
521: }
522: else {
523: ptr = x[row+mx+1];
524: }
526: dpdy = (pt - pb)/(two*hy);
527: dpdx = (pr - pl)/(two*hx);
529: /* Laplacians at each point in the 5 point stencil */
530: pblap = (pbr - two*pb + pbl)/hx2 + (p - two*pb + pbb)/hy2;
531: pllap = (p - two*pl + pll)/hx2 + (ptl - two*pl + pbl)/hy2;
532: plap = (pr - two*p + pl)/hx2 + (pt - two*p + pb)/hy2;
533: prlap = (prr - two*pr + p)/hx2 + (ptr - two*pr + pbr)/hy2;
534: ptlap = (ptr - two*pt + ptl)/hx2 + (ptt - two*pt + p)/hy2;
536: if (j > 2) {
537: val = hxhy2*(one/hy2/hy2 - rey*dpdx/hy2/(two*hy));
538: col = row - 2*mx;
539: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
540: }
541: if (i > 2) {
542: val = hxhy2*(one/hx2/hx2 + rey*dpdy/hx2/(two*hx));
543: col = row - 2;
544: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
545: }
546: if (i < mx-3) {
547: val = hxhy2*(one/hx2/hx2 - rey*dpdy/hx2/(two*hx));
548: col = row + 2;
549: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
550: }
551: if (j < my-3) {
552: val = hxhy2*(one/hy2/hy2 + rey*dpdx/hy2/(two*hy));
553: col = row + 2*mx;
554: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
555: }
556: if (i != 1 && j != 1) {
557: val = hxhy2*(two/hy2/hx2 + rey*(dpdy/hy2/(two*hx) - dpdx/hx2/(two*hy)));
558: col = row - mx - 1;
559: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
560: }
561: if (j != 1 && i != mx-2) {
562: val = hxhy2*(two/hy2/hx2 - rey*(dpdy/hy2/(two*hx) + dpdx/hx2/(two*hy)));
563: col = row - mx + 1;
564: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
565: }
566: if (j != my-2 && i != 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 != my-2 && 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 != 1) {
577: val = hxhy2*(-four*(one/hy2/hx2 + one/hy2/hy2)
578: + rey*((prlap - pllap)/(two*hx)/(two*hy)
579: + dpdx*(one/hx2 + one/hy2)/hy));
580: col = row - mx;
581: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
582: }
583: if (i != 1) {
584: val = hxhy2*(-four*(one/hy2/hx2 + one/hx2/hx2)
585: - rey*((ptlap - pblap)/(two*hx)/(two*hy)
586: + dpdy*(one/hx2 + one/hy2)/hx));
587: col = row - 1;
588: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
589: }
590: if (i != mx-2) {
591: val = hxhy2*(-four*(one/hy2/hx2 + one/hx2/hx2)
592: + rey*((ptlap - pblap)/(two*hx)/(two*hy)
593: + dpdy*(one/hx2 + one/hy2)/hx));
594: col = row + 1;
595: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
596: }
597: if (j != my-2) {
598: val = hxhy2*(-four*(one/hy2/hx2 + one/hy2/hy2)
599: - rey*((prlap - pllap)/(two*hx)/(two*hy)
600: + dpdx*(one/hx2 + one/hy2)/hy));
601: col = row + mx;
602: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
603: }
604: val = hxhy2*(two*(four/hx2/hy2 + three/hx2/hx2 + three/hy2/hy2));
605: MatSetValues(*J,1,&row,1,&row,&val,ADD_VALUES);
606: if (j == 1) {
607: val = hxhy2*(one/hy2/hy2 - rey*(dpdx/hy2/(two*hy)));
608: MatSetValues(*J,1,&row,1,&row,&val,ADD_VALUES);
609: }
610: if (i == 1) {
611: val = hxhy2*(one/hx2/hx2 + rey*(dpdy/hx2/(two*hx)));
612: MatSetValues(*J,1,&row,1,&row,&val,ADD_VALUES);
613: }
614: if (i == mx-2) {
615: val = hxhy2*(one/hx2/hx2 - rey*(dpdy/hx2/(two*hx)));
616: MatSetValues(*J,1,&row,1,&row,&val,ADD_VALUES);
617: }
618: if (j == my-2) {
619: val = hxhy2*(one/hy2/hy2 + rey*(dpdx/hy2/(two*hy)));
620: MatSetValues(*J,1,&row,1,&row,&val,ADD_VALUES);
621: }
622: }
623: }
624: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
625: VecRestoreArray(X,&x);
626: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
627: *flag = SAME_NONZERO_PATTERN;
628: return 0;
629: }