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,&its);
108: SNESGetNumberUnsuccessfulSteps(snes,&nfails);
110: PetscPrintf(PETSC_COMM_SELF,"number of Newton iterations = %d, ",its);
111: PetscPrintf(PETSC_COMM_SELF,"number of unsuccessful steps = %d\n\n",nfails);
112: VecGetArray(x,&xvalues);
113: PetscDrawTensorContour(draw,user.mx,user.my,0,0,xvalues);
114: VecRestoreArray(x,&xvalues);
116: /* Free data structures */
117: VecDestroy(x); VecDestroy(r);
118: MatDestroy(J); SNESDestroy(snes);
119: PetscDrawDestroy(draw);
120: PetscFinalize();
122: return 0;
123: }
124: /* ------------------------------------------------------------------ */
125: /* Bratu (Solid Fuel Ignition) Test Problem */
126: /* ------------------------------------------------------------------ */
128: /* -------------------- Form initial approximation ----------------- */
132: int FormInitialGuess1(AppCtx *user,Vec X)
133: {
134: int i, j, row, mx, my, ierr;
135: PetscReal lambda, temp1, temp, hx, hy, hxdhy, hydhx,sc;
136: PetscScalar *x;
138: mx = user->mx;
139: my = user->my;
140: lambda = user->param;
142: hx = 1.0 / (PetscReal)(mx-1);
143: hy = 1.0 / (PetscReal)(my-1);
144: sc = hx*hy;
145: hxdhy = hx/hy;
146: hydhx = hy/hx;
148: VecGetArray(X,&x);
149: temp1 = lambda/(lambda + 1.0);
150: for (j=0; j<my; j++) {
151: temp = (PetscReal)(PetscMin(j,my-j-1))*hy;
152: for (i=0; i<mx; i++) {
153: row = i + j*mx;
154: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
155: x[row] = 0.0;
156: continue;
157: }
158: x[row] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));
159: }
160: }
161: VecRestoreArray(X,&x);
162: return 0;
163: }
164: /* -------------------- Evaluate Function F(x) --------------------- */
165:
168: int FormFunction1(SNES snes,Vec X,Vec F,void *ptr)
169: {
170: AppCtx *user = (AppCtx*)ptr;
171: int ierr, i, j, row, mx, my;
172: PetscReal two = 2.0, one = 1.0, lambda,hx, hy, hxdhy, hydhx;
173: PetscScalar ut, ub, ul, ur, u, uxx, uyy, sc,*x,*f;
175: mx = user->mx;
176: my = user->my;
177: lambda = user->param;
179: hx = one / (PetscReal)(mx-1);
180: hy = one / (PetscReal)(my-1);
181: sc = hx*hy;
182: hxdhy = hx/hy;
183: hydhx = hy/hx;
185: VecGetArray(X,&x);
186: VecGetArray(F,&f);
187: for (j=0; j<my; j++) {
188: for (i=0; i<mx; i++) {
189: row = i + j*mx;
190: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
191: f[row] = x[row];
192: continue;
193: }
194: u = x[row];
195: ub = x[row - mx];
196: ul = x[row - 1];
197: ut = x[row + mx];
198: ur = x[row + 1];
199: uxx = (-ur + two*u - ul)*hydhx;
200: uyy = (-ut + two*u - ub)*hxdhy;
201: f[row] = uxx + uyy - sc*lambda*exp(u);
202: }
203: }
204: VecRestoreArray(X,&x);
205: VecRestoreArray(F,&f);
206: return 0;
207: }
208: /* -------------------- Evaluate Jacobian F'(x) -------------------- */
212: int FormJacobian1(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
213: {
214: AppCtx *user = (AppCtx*)ptr;
215: Mat jac = *J;
216: int i, j, row, mx, my, col[5], ierr;
217: PetscScalar two = 2.0, one = 1.0, lambda, v[5],sc, *x;
218: PetscReal hx, hy, hxdhy, hydhx;
221: mx = user->mx;
222: my = user->my;
223: lambda = user->param;
225: hx = 1.0 / (PetscReal)(mx-1);
226: hy = 1.0 / (PetscReal)(my-1);
227: sc = hx*hy;
228: hxdhy = hx/hy;
229: hydhx = hy/hx;
231: VecGetArray(X,&x);
232: for (j=0; j<my; j++) {
233: for (i=0; i<mx; i++) {
234: row = i + j*mx;
235: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
236: MatSetValues(jac,1,&row,1,&row,&one,INSERT_VALUES);
237: continue;
238: }
239: v[0] = -hxdhy; col[0] = row - mx;
240: v[1] = -hydhx; col[1] = row - 1;
241: v[2] = two*(hydhx + hxdhy) - sc*lambda*exp(x[row]); col[2] = row;
242: v[3] = -hydhx; col[3] = row + 1;
243: v[4] = -hxdhy; col[4] = row + mx;
244: MatSetValues(jac,1,&row,5,col,v,INSERT_VALUES);
245: }
246: }
247: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
248: VecRestoreArray(X,&x);
249: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
250: *flag = SAME_NONZERO_PATTERN;
251: return 0;
252: }
253: /* ------------------------------------------------------------------ */
254: /* Driven Cavity Test Problem */
255: /* ------------------------------------------------------------------ */
257: /* -------------------- Form initial approximation ----------------- */
261: int FormInitialGuess2(AppCtx *user,Vec X)
262: {
263: int ierr, i, j, row, mx, my;
264: PetscScalar xx,yy,*x;
265: PetscReal hx, hy;
267: mx = user->mx;
268: my = user->my;
270: /* Test for invalid input parameters */
271: if ((mx <= 0) || (my <= 0)) SETERRQ(1,0);
273: hx = 1.0 / (PetscReal)(mx-1);
274: hy = 1.0 / (PetscReal)(my-1);
276: VecGetArray(X,&x);
277: yy = 0.0;
278: for (j=0; j<my; j++) {
279: xx = 0.0;
280: for (i=0; i<mx; i++) {
281: row = i + j*mx;
282: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
283: x[row] = 0.0;
284: }
285: else {
286: x[row] = - xx*(1.0 - xx)*yy*(1.0 - yy);
287: }
288: xx = xx + hx;
289: }
290: yy = yy + hy;
291: }
292: VecRestoreArray(X,&x);
293: return 0;
294: }
295: /* -------------------- Evaluate Function F(x) --------------------- */
299: int FormFunction2(SNES snes,Vec X,Vec F,void *pptr)
300: {
301: AppCtx *user = (AppCtx*)pptr;
302: int i, j, row, mx, my, ierr;
303: PetscScalar two = 2.0, zero = 0.0, pb, pbb,pbr, pl,pll,p,pr,prr;
304: PetscScalar ptl,pt,ptt,dpdy,dpdx,pblap,ptlap,rey,pbl,ptr,pllap,plap,prlap;
305: PetscScalar *x,*f, hx2, hy2, hxhy2;
306: PetscReal hx, hy;
308: mx = user->mx;
309: my = user->my;
310: hx = 1.0 / (PetscReal)(mx-1);
311: hy = 1.0 / (PetscReal)(my-1);
312: hx2 = hx*hx;
313: hy2 = hy*hy;
314: hxhy2 = hx2*hy2;
315: rey = user->param;
317: VecGetArray(X,&x);
318: VecGetArray(F,&f);
319: for (j=0; j<my; j++) {
320: for (i=0; i<mx; i++) {
321: row = i + j*mx;
322: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
323: f[row] = x[row];
324: continue;
325: }
326: if (i == 1 || j == 1) {
327: pbl = zero;
328: }
329: else {
330: pbl = x[row-mx-1];
331: }
332: if (j == 1) {
333: pb = zero;
334: pbb = x[row];
335: }
336: else if (j == 2) {
337: pb = x[row-mx];
338: pbb = zero;
339: }
340: else {
341: pb = x[row-mx];
342: pbb = x[row-2*mx];
343: }
344: if (j == 1 || i == mx-2) {
345: pbr = zero;
346: }
347: else {
348: pbr = x[row-mx+1];
349: }
350: if (i == 1) {
351: pl = zero;
352: pll = x[row];
353: }
354: else if (i == 2) {
355: pl = x[row-1];
356: pll = zero;
357: }
358: else {
359: pl = x[row-1];
360: pll = x[row-2];
361: }
362: p = x[row];
363: if (i == mx-3) {
364: pr = x[row+1];
365: prr = zero;
366: }
367: else if (i == mx-2) {
368: pr = zero;
369: prr = x[row];
370: }
371: else {
372: pr = x[row+1];
373: prr = x[row+2];
374: }
375: if (j == my-2 || i == 1) {
376: ptl = zero;
377: }
378: else {
379: ptl = x[row+mx-1];
380: }
381: if (j == my-3) {
382: pt = x[row+mx];
383: ptt = zero;
384: }
385: else if (j == my-2) {
386: pt = zero;
387: ptt = x[row] + two*hy;
388: }
389: else {
390: pt = x[row+mx];
391: ptt = x[row+2*mx];
392: }
393: if (j == my-2 || i == mx-2) {
394: ptr = zero;
395: }
396: else {
397: ptr = x[row+mx+1];
398: }
400: dpdy = (pt - pb)/(two*hy);
401: dpdx = (pr - pl)/(two*hx);
403: /* Laplacians at each point in the 5 point stencil */
404: pblap = (pbr - two*pb + pbl)/hx2 + (p - two*pb + pbb)/hy2;
405: pllap = (p - two*pl + pll)/hx2 + (ptl - two*pl + pbl)/hy2;
406: plap = (pr - two*p + pl)/hx2 + (pt - two*p + pb)/hy2;
407: prlap = (prr - two*pr + p)/hx2 + (ptr - two*pr + pbr)/hy2;
408: ptlap = (ptr - two*pt + ptl)/hx2 + (ptt - two*pt + p)/hy2;
410: f[row] = hxhy2*((prlap - two*plap + pllap)/hx2
411: + (ptlap - two*plap + pblap)/hy2
412: - rey*(dpdy*(prlap - pllap)/(two*hx)
413: - dpdx*(ptlap - pblap)/(two*hy)));
414: }
415: }
416: VecRestoreArray(X,&x);
417: VecRestoreArray(F,&f);
418: return 0;
419: }
420: /* -------------------- Evaluate Jacobian F'(x) -------------------- */
424: int FormJacobian2(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *pptr)
425: {
426: AppCtx *user = (AppCtx*)pptr;
427: int i, j, row, mx, my, col, ierr;
428: PetscScalar two = 2.0, one = 1.0, zero = 0.0, pb, pbb,pbr, pl,pll,p,pr,prr;
429: PetscScalar ptl,pt,ptt,dpdy,dpdx,pblap,ptlap,rey,pbl,ptr,pllap,plap,prlap;
430: PetscScalar val,four = 4.0, three = 3.0,*x;
431: PetscReal hx, hy,hx2, hy2, hxhy2;
433: mx = user->mx;
434: my = user->my;
435: hx = 1.0 / (PetscReal)(mx-1);
436: hy = 1.0 / (PetscReal)(my-1);
437: hx2 = hx*hx;
438: hy2 = hy*hy;
439: hxhy2 = hx2*hy2;
440: rey = user->param;
442: MatZeroEntries(*J);
443: VecGetArray(X,&x);
444: for (j=0; j<my; j++) {
445: for (i=0; i<mx; i++) {
446: row = i + j*mx;
447: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
448: MatSetValues(*J,1,&row,1,&row,&one,ADD_VALUES);
449: continue;
450: }
451: if (i == 1 || j == 1) {
452: pbl = zero;
453: }
454: else {
455: pbl = x[row-mx-1];
456: }
457: if (j == 1) {
458: pb = zero;
459: pbb = x[row];
460: }
461: else if (j == 2) {
462: pb = x[row-mx];
463: pbb = zero;
464: }
465: else {
466: pb = x[row-mx];
467: pbb = x[row-2*mx];
468: }
469: if (j == 1 || i == mx-2) {
470: pbr = zero;
471: }
472: else {
473: pbr = x[row-mx+1];
474: }
475: if (i == 1) {
476: pl = zero;
477: pll = x[row];
478: }
479: else if (i == 2) {
480: pl = x[row-1];
481: pll = zero;
482: }
483: else {
484: pl = x[row-1];
485: pll = x[row-2];
486: }
487: p = x[row];
488: if (i == mx-3) {
489: pr = x[row+1];
490: prr = zero;
491: }
492: else if (i == mx-2) {
493: pr = zero;
494: prr = x[row];
495: }
496: else {
497: pr = x[row+1];
498: prr = x[row+2];
499: }
500: if (j == my-2 || i == 1) {
501: ptl = zero;
502: }
503: else {
504: ptl = x[row+mx-1];
505: }
506: if (j == my-3) {
507: pt = x[row+mx];
508: ptt = zero;
509: }
510: else if (j == my-2) {
511: pt = zero;
512: ptt = x[row] + two*hy;
513: }
514: else {
515: pt = x[row+mx];
516: ptt = x[row+2*mx];
517: }
518: if (j == my-2 || i == mx-2) {
519: ptr = zero;
520: }
521: else {
522: ptr = x[row+mx+1];
523: }
525: dpdy = (pt - pb)/(two*hy);
526: dpdx = (pr - pl)/(two*hx);
528: /* Laplacians at each point in the 5 point stencil */
529: pblap = (pbr - two*pb + pbl)/hx2 + (p - two*pb + pbb)/hy2;
530: pllap = (p - two*pl + pll)/hx2 + (ptl - two*pl + pbl)/hy2;
531: plap = (pr - two*p + pl)/hx2 + (pt - two*p + pb)/hy2;
532: prlap = (prr - two*pr + p)/hx2 + (ptr - two*pr + pbr)/hy2;
533: ptlap = (ptr - two*pt + ptl)/hx2 + (ptt - two*pt + p)/hy2;
535: if (j > 2) {
536: val = hxhy2*(one/hy2/hy2 - rey*dpdx/hy2/(two*hy));
537: col = row - 2*mx;
538: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
539: }
540: if (i > 2) {
541: val = hxhy2*(one/hx2/hx2 + rey*dpdy/hx2/(two*hx));
542: col = row - 2;
543: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
544: }
545: if (i < mx-3) {
546: val = hxhy2*(one/hx2/hx2 - rey*dpdy/hx2/(two*hx));
547: col = row + 2;
548: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
549: }
550: if (j < my-3) {
551: val = hxhy2*(one/hy2/hy2 + rey*dpdx/hy2/(two*hy));
552: col = row + 2*mx;
553: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
554: }
555: if (i != 1 && j != 1) {
556: val = hxhy2*(two/hy2/hx2 + rey*(dpdy/hy2/(two*hx) - dpdx/hx2/(two*hy)));
557: col = row - mx - 1;
558: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
559: }
560: if (j != 1 && i != mx-2) {
561: val = hxhy2*(two/hy2/hx2 - rey*(dpdy/hy2/(two*hx) + dpdx/hx2/(two*hy)));
562: col = row - mx + 1;
563: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
564: }
565: if (j != my-2 && i != 1) {
566: val = hxhy2*(two/hy2/hx2 + rey*(dpdy/hy2/(two*hx) + dpdx/hx2/(two*hy)));
567: col = row + mx - 1;
568: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
569: }
570: if (j != my-2 && i != mx-2) {
571: val = hxhy2*(two/hy2/hx2 - rey*(dpdy/hy2/(two*hx) - dpdx/hx2/(two*hy)));
572: col = row + mx + 1;
573: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
574: }
575: if (j != 1) {
576: val = hxhy2*(-four*(one/hy2/hx2 + one/hy2/hy2)
577: + rey*((prlap - pllap)/(two*hx)/(two*hy)
578: + dpdx*(one/hx2 + one/hy2)/hy));
579: col = row - mx;
580: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
581: }
582: if (i != 1) {
583: val = hxhy2*(-four*(one/hy2/hx2 + one/hx2/hx2)
584: - rey*((ptlap - pblap)/(two*hx)/(two*hy)
585: + dpdy*(one/hx2 + one/hy2)/hx));
586: col = row - 1;
587: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
588: }
589: if (i != mx-2) {
590: val = hxhy2*(-four*(one/hy2/hx2 + one/hx2/hx2)
591: + rey*((ptlap - pblap)/(two*hx)/(two*hy)
592: + dpdy*(one/hx2 + one/hy2)/hx));
593: col = row + 1;
594: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
595: }
596: if (j != my-2) {
597: val = hxhy2*(-four*(one/hy2/hx2 + one/hy2/hy2)
598: - rey*((prlap - pllap)/(two*hx)/(two*hy)
599: + dpdx*(one/hx2 + one/hy2)/hy));
600: col = row + mx;
601: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
602: }
603: val = hxhy2*(two*(four/hx2/hy2 + three/hx2/hx2 + three/hy2/hy2));
604: MatSetValues(*J,1,&row,1,&row,&val,ADD_VALUES);
605: if (j == 1) {
606: val = hxhy2*(one/hy2/hy2 - rey*(dpdx/hy2/(two*hy)));
607: MatSetValues(*J,1,&row,1,&row,&val,ADD_VALUES);
608: }
609: if (i == 1) {
610: val = hxhy2*(one/hx2/hx2 + rey*(dpdy/hx2/(two*hx)));
611: MatSetValues(*J,1,&row,1,&row,&val,ADD_VALUES);
612: }
613: if (i == mx-2) {
614: val = hxhy2*(one/hx2/hx2 - rey*(dpdy/hx2/(two*hx)));
615: MatSetValues(*J,1,&row,1,&row,&val,ADD_VALUES);
616: }
617: if (j == my-2) {
618: val = hxhy2*(one/hy2/hy2 + rey*(dpdx/hy2/(two*hy)));
619: MatSetValues(*J,1,&row,1,&row,&val,ADD_VALUES);
620: }
621: }
622: }
623: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
624: VecRestoreArray(X,&x);
625: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
626: *flag = SAME_NONZERO_PATTERN;
627: return 0;
628: }