Actual source code: ex17.c
petsc-3.8.4 2018-03-24
1: static const char help[] = "Newton's method to solve a two-variable system, sequentially.\n"
2: "The same problem is solved twice - i) fully assembled system + ii) block system\n\n";
4: /*T
5: Concepts: SNES^basic uniprocessor example, block objects
6: Processors: 1
7: T*/
9: /*
10: Include "petscsnes.h" so that we can use SNES solvers. Note that this
11: file automatically includes:
12: petscsys.h - base PETSc routines petscvec.h - vectors
13: petscsys.h - system routines petscmat.h - matrices
14: petscis.h - index sets petscksp.h - Krylov subspace methods
15: petscviewer.h - viewers petscpc.h - preconditioners
16: petscksp.h - linear solvers
17: */
18: #include <petscsnes.h>
20: /*
21: This example is block version of the test found at
22: ${PETSC_DIR}/src/snes/examples/tutorials/ex1.c
23: In this test we replace the Jacobian systems
24: [J]{x} = {F}
25: where
27: [J] = (j_00, j_01), {x} = (x_0, x_1)^T, {F} = (f_0, f_1)^T
28: (j_10, j_11)
29: where [J] \in \mathbb^{2 \times 2}, {x},{F} \in \mathbb^{2 \times 1},
31: with a block system in which each block is of length 1.
32: i.e. The block system is thus
34: [J] = ([j00], [j01]), {x} = ({x0}, {x1})^T, {F} = ({f0}, {f1})^T
35: ([j10], [j11])
36: where
37: [j00], [j01], [j10], [j11] \in \mathbb^{1 \times 1}
38: {x0}, {x1}, {f0}, {f1} \in \mathbb^{1 \times 1}
40: In practice we would not bother defing blocks of size one, and would instead assemble the
41: full system. This is just a simple test to illustrate how to manipulate the blocks and
42: to confirm the implementation is correct.
43: */
45: /*
46: User-defined routines
47: */
48: static PetscErrorCode FormJacobian1(SNES,Vec,Mat,Mat,void*);
49: static PetscErrorCode FormFunction1(SNES,Vec,Vec,void*);
50: static PetscErrorCode FormJacobian2(SNES,Vec,Mat,Mat,void*);
51: static PetscErrorCode FormFunction2(SNES,Vec,Vec,void*);
52: static PetscErrorCode FormJacobian1_block(SNES,Vec,Mat,Mat,void*);
53: static PetscErrorCode FormFunction1_block(SNES,Vec,Vec,void*);
54: static PetscErrorCode FormJacobian2_block(SNES,Vec,Mat,Mat,void*);
55: static PetscErrorCode FormFunction2_block(SNES,Vec,Vec,void*);
58: static PetscErrorCode assembled_system(void)
59: {
60: SNES snes; /* nonlinear solver context */
61: KSP ksp; /* linear solver context */
62: PC pc; /* preconditioner context */
63: Vec x,r; /* solution, residual vectors */
64: Mat J; /* Jacobian matrix */
66: PetscInt its;
67: PetscScalar pfive = .5,*xx;
68: PetscBool flg;
71: PetscPrintf(PETSC_COMM_WORLD, "\n\n========================= Assembled system =========================\n\n");
73: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
74: Create nonlinear solver context
75: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
77: SNESCreate(PETSC_COMM_WORLD,&snes);
79: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
80: Create matrix and vector data structures; set corresponding routines
81: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
83: /*
84: Create vectors for solution and nonlinear function
85: */
86: VecCreateSeq(PETSC_COMM_SELF,2,&x);
87: VecDuplicate(x,&r);
89: /*
90: Create Jacobian matrix data structure
91: */
92: MatCreate(PETSC_COMM_SELF,&J);
93: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);
94: MatSetFromOptions(J);
95: MatSetUp(J);
97: PetscOptionsHasName(NULL,NULL,"-hard",&flg);
98: if (!flg) {
99: /*
100: Set function evaluation routine and vector.
101: */
102: SNESSetFunction(snes,r,FormFunction1,NULL);
104: /*
105: Set Jacobian matrix data structure and Jacobian evaluation routine
106: */
107: SNESSetJacobian(snes,J,J,FormJacobian1,NULL);
108: } else {
109: SNESSetFunction(snes,r,FormFunction2,NULL);
110: SNESSetJacobian(snes,J,J,FormJacobian2,NULL);
111: }
113: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: Customize nonlinear solver; set runtime options
115: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
117: /*
118: Set linear solver defaults for this problem. By extracting the
119: KSP, KSP, and PC contexts from the SNES context, we can then
120: directly call any KSP, KSP, and PC routines to set various options.
121: */
122: SNESGetKSP(snes,&ksp);
123: KSPGetPC(ksp,&pc);
124: PCSetType(pc,PCNONE);
125: KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);
127: /*
128: Set SNES/KSP/KSP/PC runtime options, e.g.,
129: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
130: These options will override those specified above as long as
131: SNESSetFromOptions() is called _after_ any other customization
132: routines.
133: */
134: SNESSetFromOptions(snes);
136: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137: Evaluate initial guess; then solve nonlinear system
138: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
139: if (!flg) {
140: VecSet(x,pfive);
141: } else {
142: VecGetArray(x,&xx);
143: xx[0] = 2.0; xx[1] = 3.0;
144: VecRestoreArray(x,&xx);
145: }
146: /*
147: Note: The user should initialize the vector, x, with the initial guess
148: for the nonlinear solver prior to calling SNESSolve(). In particular,
149: to employ an initial guess of zero, the user should explicitly set
150: this vector to zero by calling VecSet().
151: */
153: SNESSolve(snes,NULL,x);
154: SNESGetIterationNumber(snes,&its);
155: if (flg) {
156: Vec f;
157: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
158: SNESGetFunction(snes,&f,0,0);
159: VecView(r,PETSC_VIEWER_STDOUT_WORLD);
160: }
162: PetscPrintf(PETSC_COMM_SELF,"number of SNES iterations = %D\n\n",its);
164: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
165: Free work space. All PETSc objects should be destroyed when they
166: are no longer needed.
167: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
169: VecDestroy(&x); VecDestroy(&r);
170: MatDestroy(&J); SNESDestroy(&snes);
171: return(0);
172: }
174: /* ------------------------------------------------------------------- */
175: /*
176: FormFunction1 - Evaluates nonlinear function, F(x).
178: Input Parameters:
179: . snes - the SNES context
180: . x - input vector
181: . dummy - optional user-defined context (not used here)
183: Output Parameter:
184: . f - function vector
185: */
186: static PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *dummy)
187: {
188: PetscErrorCode ierr;
189: const PetscScalar *xx;
190: PetscScalar *ff;
193: /*
194: Get pointers to vector data.
195: - For default PETSc vectors, VecGetArray() returns a pointer to
196: the data array. Otherwise, the routine is implementation dependent.
197: - You MUST call VecRestoreArray() when you no longer need access to
198: the array.
199: */
200: VecGetArrayRead(x,&xx);
201: VecGetArray(f,&ff);
203: /*
204: Compute function
205: */
206: ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0;
207: ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0;
210: /*
211: Restore vectors
212: */
213: VecRestoreArrayRead(x,&xx);
214: VecRestoreArray(f,&ff);
215: return(0);
216: }
217: /* ------------------------------------------------------------------- */
218: /*
219: FormJacobian1 - Evaluates Jacobian matrix.
221: Input Parameters:
222: . snes - the SNES context
223: . x - input vector
224: . dummy - optional user-defined context (not used here)
226: Output Parameters:
227: . jac - Jacobian matrix
228: . B - optionally different preconditioning matrix
229: . flag - flag indicating matrix structure
230: */
231: static PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
232: {
233: const PetscScalar *xx;
234: PetscScalar A[4];
235: PetscErrorCode ierr;
236: PetscInt idx[2] = {0,1};
239: /*
240: Get pointer to vector data
241: */
242: VecGetArrayRead(x,&xx);
244: /*
245: Compute Jacobian entries and insert into matrix.
246: - Since this is such a small problem, we set all entries for
247: the matrix at once.
248: */
249: A[0] = 2.0*xx[0] + xx[1]; A[1] = xx[0];
250: A[2] = xx[1]; A[3] = xx[0] + 2.0*xx[1];
251: MatSetValues(jac,2,idx,2,idx,A,INSERT_VALUES);
253: /*
254: Restore vector
255: */
256: VecRestoreArrayRead(x,&xx);
258: /*
259: Assemble matrix
260: */
261: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
262: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
263: return(0);
264: }
267: /* ------------------------------------------------------------------- */
268: static PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
269: {
270: PetscErrorCode ierr;
271: const PetscScalar *xx;
272: PetscScalar *ff;
275: /*
276: Get pointers to vector data.
277: - For default PETSc vectors, VecGetArray() returns a pointer to
278: the data array. Otherwise, the routine is implementation dependent.
279: - You MUST call VecRestoreArray() when you no longer need access to
280: the array.
281: */
282: VecGetArrayRead(x,&xx);
283: VecGetArray(f,&ff);
285: /*
286: Compute function
287: */
288: ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
289: ff[1] = xx[1];
291: /*
292: Restore vectors
293: */
294: VecRestoreArrayRead(x,&xx);
295: VecRestoreArray(f,&ff);
296: return(0);
297: }
299: /* ------------------------------------------------------------------- */
300: static PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
301: {
302: const PetscScalar *xx;
303: PetscScalar A[4];
304: PetscErrorCode ierr;
305: PetscInt idx[2] = {0,1};
308: /*
309: Get pointer to vector data
310: */
311: VecGetArrayRead(x,&xx);
313: /*
314: Compute Jacobian entries and insert into matrix.
315: - Since this is such a small problem, we set all entries for
316: the matrix at once.
317: */
318: A[0] = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
319: A[2] = 0.0; A[3] = 1.0;
320: MatSetValues(jac,2,idx,2,idx,A,INSERT_VALUES);
322: /*
323: Restore vector
324: */
325: VecRestoreArrayRead(x,&xx);
327: /*
328: Assemble matrix
329: */
330: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
331: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
332: return(0);
333: }
335: static int block_system(void)
336: {
337: SNES snes; /* nonlinear solver context */
338: KSP ksp; /* linear solver context */
339: PC pc; /* preconditioner context */
340: Vec x,r; /* solution, residual vectors */
341: Mat J; /* Jacobian matrix */
343: PetscInt its;
344: PetscScalar pfive = .5;
345: PetscBool flg;
347: Mat j11, j12, j21, j22;
348: Vec x1, x2, r1, r2;
349: Vec bv;
350: Vec bx[2];
351: Mat bA[2][2];
354: PetscPrintf(PETSC_COMM_WORLD, "\n\n========================= Block system =========================\n\n");
356: SNESCreate(PETSC_COMM_WORLD,&snes);
358: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
359: Create matrix and vector data structures; set corresponding routines
360: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
362: /*
363: Create sub vectors for solution and nonlinear function
364: */
365: VecCreateSeq(PETSC_COMM_SELF,1,&x1);
366: VecDuplicate(x1,&r1);
368: VecCreateSeq(PETSC_COMM_SELF,1,&x2);
369: VecDuplicate(x2,&r2);
371: /*
372: Create the block vectors
373: */
374: bx[0] = x1;
375: bx[1] = x2;
376: VecCreateNest(PETSC_COMM_WORLD,2,NULL,bx,&x);
377: VecAssemblyBegin(x);
378: VecAssemblyEnd(x);
379: VecDestroy(&x1);
380: VecDestroy(&x2);
382: bx[0] = r1;
383: bx[1] = r2;
384: VecCreateNest(PETSC_COMM_WORLD,2,NULL,bx,&r);
385: VecDestroy(&r1);
386: VecDestroy(&r2);
387: VecAssemblyBegin(r);
388: VecAssemblyEnd(r);
390: /*
391: Create sub Jacobian matrix data structure
392: */
393: MatCreate(PETSC_COMM_WORLD, &j11);
394: MatSetSizes(j11, 1, 1, 1, 1);
395: MatSetType(j11, MATSEQAIJ);
396: MatSetUp(j11);
398: MatCreate(PETSC_COMM_WORLD, &j12);
399: MatSetSizes(j12, 1, 1, 1, 1);
400: MatSetType(j12, MATSEQAIJ);
401: MatSetUp(j12);
403: MatCreate(PETSC_COMM_WORLD, &j21);
404: MatSetSizes(j21, 1, 1, 1, 1);
405: MatSetType(j21, MATSEQAIJ);
406: MatSetUp(j21);
408: MatCreate(PETSC_COMM_WORLD, &j22);
409: MatSetSizes(j22, PETSC_DECIDE, PETSC_DECIDE, 1, 1);
410: MatSetType(j22, MATSEQAIJ);
411: MatSetUp(j22);
412: /*
413: Create block Jacobian matrix data structure
414: */
415: bA[0][0] = j11;
416: bA[0][1] = j12;
417: bA[1][0] = j21;
418: bA[1][1] = j22;
420: MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,&bA[0][0],&J);
421: MatSetUp(J);
422: MatNestSetVecType(J,VECNEST);
423: MatDestroy(&j11);
424: MatDestroy(&j12);
425: MatDestroy(&j21);
426: MatDestroy(&j22);
428: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
429: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
431: PetscOptionsHasName(NULL,NULL,"-hard",&flg);
432: if (!flg) {
433: /*
434: Set function evaluation routine and vector.
435: */
436: SNESSetFunction(snes,r,FormFunction1_block,NULL);
438: /*
439: Set Jacobian matrix data structure and Jacobian evaluation routine
440: */
441: SNESSetJacobian(snes,J,J,FormJacobian1_block,NULL);
442: } else {
443: SNESSetFunction(snes,r,FormFunction2_block,NULL);
444: SNESSetJacobian(snes,J,J,FormJacobian2_block,NULL);
445: }
447: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
448: Customize nonlinear solver; set runtime options
449: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
451: /*
452: Set linear solver defaults for this problem. By extracting the
453: KSP, KSP, and PC contexts from the SNES context, we can then
454: directly call any KSP, KSP, and PC routines to set various options.
455: */
456: SNESGetKSP(snes,&ksp);
457: KSPGetPC(ksp,&pc);
458: PCSetType(pc,PCNONE);
459: KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);
461: /*
462: Set SNES/KSP/KSP/PC runtime options, e.g.,
463: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
464: These options will override those specified above as long as
465: SNESSetFromOptions() is called _after_ any other customization
466: routines.
467: */
468: SNESSetFromOptions(snes);
470: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
471: Evaluate initial guess; then solve nonlinear system
472: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
473: if (!flg) {
474: VecSet(x,pfive);
475: } else {
476: Vec *vecs;
477: VecNestGetSubVecs(x, NULL, &vecs);
478: bv = vecs[0];
479: /* VecBlockGetSubVec(x, 0, &bv); */
480: VecSetValue(bv, 0, 2.0, INSERT_VALUES); /* xx[0] = 2.0; */
481: VecAssemblyBegin(bv);
482: VecAssemblyEnd(bv);
484: /* VecBlockGetSubVec(x, 1, &bv); */
485: bv = vecs[1];
486: VecSetValue(bv, 0, 3.0, INSERT_VALUES); /* xx[1] = 3.0; */
487: VecAssemblyBegin(bv);
488: VecAssemblyEnd(bv);
489: }
490: /*
491: Note: The user should initialize the vector, x, with the initial guess
492: for the nonlinear solver prior to calling SNESSolve(). In particular,
493: to employ an initial guess of zero, the user should explicitly set
494: this vector to zero by calling VecSet().
495: */
496: SNESSolve(snes,NULL,x);
497: SNESGetIterationNumber(snes,&its);
498: if (flg) {
499: Vec f;
500: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
501: SNESGetFunction(snes,&f,0,0);
502: VecView(r,PETSC_VIEWER_STDOUT_WORLD);
503: }
505: PetscPrintf(PETSC_COMM_SELF,"number of SNES iterations = %D\n\n",its);
507: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
508: Free work space. All PETSc objects should be destroyed when they
509: are no longer needed.
510: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
511: VecDestroy(&x); VecDestroy(&r);
512: MatDestroy(&J); SNESDestroy(&snes);
513: return(0);
514: }
516: /* ------------------------------------------------------------------- */
517: static PetscErrorCode FormFunction1_block(SNES snes,Vec x,Vec f,void *dummy)
518: {
519: Vec *xx, *ff, x1,x2, f1,f2;
520: PetscScalar ff_0, ff_1;
521: PetscScalar xx_0, xx_1;
522: PetscInt index,nb;
526: /* get blocks for function */
527: VecNestGetSubVecs(f, &nb, &ff);
528: f1 = ff[0]; f2 = ff[1];
530: /* get blocks for solution */
531: VecNestGetSubVecs(x, &nb, &xx);
532: x1 = xx[0]; x2 = xx[1];
534: /* get solution values */
535: index = 0;
536: VecGetValues(x1,1, &index, &xx_0);
537: VecGetValues(x2,1, &index, &xx_1);
539: /* Compute function */
540: ff_0 = xx_0*xx_0 + xx_0*xx_1 - 3.0;
541: ff_1 = xx_0*xx_1 + xx_1*xx_1 - 6.0;
543: /* set function values */
544: VecSetValue(f1, index, ff_0, INSERT_VALUES);
546: VecSetValue(f2, index, ff_1, INSERT_VALUES);
548: VecAssemblyBegin(f);
549: VecAssemblyEnd(f);
550: return(0);
551: }
553: /* ------------------------------------------------------------------- */
554: static PetscErrorCode FormJacobian1_block(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
555: {
556: Vec *xx, x1,x2;
557: PetscScalar xx_0, xx_1;
558: PetscInt index,nb;
559: PetscScalar A_00, A_01, A_10, A_11;
560: Mat j11, j12, j21, j22;
561: Mat **mats;
565: /* get blocks for solution */
566: VecNestGetSubVecs(x, &nb, &xx);
567: x1 = xx[0]; x2 = xx[1];
569: /* get solution values */
570: index = 0;
571: VecGetValues(x1,1, &index, &xx_0);
572: VecGetValues(x2,1, &index, &xx_1);
574: /* get block matrices */
575: MatNestGetSubMats(jac,NULL,NULL,&mats);
576: j11 = mats[0][0];
577: j12 = mats[0][1];
578: j21 = mats[1][0];
579: j22 = mats[1][1];
581: /* compute jacobian entries */
582: A_00 = 2.0*xx_0 + xx_1;
583: A_01 = xx_0;
584: A_10 = xx_1;
585: A_11 = xx_0 + 2.0*xx_1;
587: /* set jacobian values */
588: MatSetValue(j11, 0,0, A_00, INSERT_VALUES);
589: MatSetValue(j12, 0,0, A_01, INSERT_VALUES);
590: MatSetValue(j21, 0,0, A_10, INSERT_VALUES);
591: MatSetValue(j22, 0,0, A_11, INSERT_VALUES);
593: /* Assemble sub matrix */
594: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
595: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
596: return(0);
597: }
599: /* ------------------------------------------------------------------- */
600: static PetscErrorCode FormFunction2_block(SNES snes,Vec x,Vec f,void *dummy)
601: {
602: PetscErrorCode ierr;
603: PetscScalar *ff;
604: const PetscScalar *xx;
607: /*
608: Get pointers to vector data.
609: - For default PETSc vectors, VecGetArray() returns a pointer to
610: the data array. Otherwise, the routine is implementation dependent.
611: - You MUST call VecRestoreArray() when you no longer need access to
612: the array.
613: */
614: VecGetArrayRead(x,&xx);
615: VecGetArray(f,&ff);
617: /*
618: Compute function
619: */
620: ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
621: ff[1] = xx[1];
623: /*
624: Restore vectors
625: */
626: VecRestoreArrayRead(x,&xx);
627: VecRestoreArray(f,&ff);
628: return(0);
629: }
631: /* ------------------------------------------------------------------- */
632: static PetscErrorCode FormJacobian2_block(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
633: {
634: const PetscScalar *xx;
635: PetscScalar A[4];
636: PetscErrorCode ierr;
637: PetscInt idx[2] = {0,1};
640: /*
641: Get pointer to vector data
642: */
643: VecGetArrayRead(x,&xx);
645: /*
646: Compute Jacobian entries and insert into matrix.
647: - Since this is such a small problem, we set all entries for
648: the matrix at once.
649: */
650: A[0] = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
651: A[2] = 0.0; A[3] = 1.0;
652: MatSetValues(jac,2,idx,2,idx,A,INSERT_VALUES);
654: /*
655: Restore vector
656: */
657: VecRestoreArrayRead(x,&xx);
659: /*
660: Assemble matrix
661: */
662: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
663: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
664: return(0);
665: }
667: int main(int argc,char **argv)
668: {
669: PetscMPIInt size;
672: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
674: MPI_Comm_size(PETSC_COMM_WORLD,&size);
675: if (size != 1) SETERRQ(PETSC_COMM_WORLD, 1,"This is a uniprocessor example only!");
677: assembled_system();
679: block_system();
681: PetscFinalize();
682: return 0;
683: }