Actual source code: ex17.c
petsc-3.6.4 2016-04-12
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*);
60: static PetscErrorCode assembled_system(void)
61: {
62: SNES snes; /* nonlinear solver context */
63: KSP ksp; /* linear solver context */
64: PC pc; /* preconditioner context */
65: Vec x,r; /* solution, residual vectors */
66: Mat J; /* Jacobian matrix */
68: PetscInt its;
69: PetscScalar pfive = .5,*xx;
70: PetscBool flg;
73: PetscPrintf(PETSC_COMM_WORLD, "\n\n========================= Assembled system =========================\n\n");
75: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
76: Create nonlinear solver context
77: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
79: SNESCreate(PETSC_COMM_WORLD,&snes);
81: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
82: Create matrix and vector data structures; set corresponding routines
83: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
85: /*
86: Create vectors for solution and nonlinear function
87: */
88: VecCreateSeq(PETSC_COMM_SELF,2,&x);
89: VecDuplicate(x,&r);
91: /*
92: Create Jacobian matrix data structure
93: */
94: MatCreate(PETSC_COMM_SELF,&J);
95: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);
96: MatSetFromOptions(J);
97: MatSetUp(J);
99: PetscOptionsHasName(NULL,"-hard",&flg);
100: if (!flg) {
101: /*
102: Set function evaluation routine and vector.
103: */
104: SNESSetFunction(snes,r,FormFunction1,NULL);
106: /*
107: Set Jacobian matrix data structure and Jacobian evaluation routine
108: */
109: SNESSetJacobian(snes,J,J,FormJacobian1,NULL);
110: } else {
111: SNESSetFunction(snes,r,FormFunction2,NULL);
112: SNESSetJacobian(snes,J,J,FormJacobian2,NULL);
113: }
115: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: Customize nonlinear solver; set runtime options
117: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
119: /*
120: Set linear solver defaults for this problem. By extracting the
121: KSP, KSP, and PC contexts from the SNES context, we can then
122: directly call any KSP, KSP, and PC routines to set various options.
123: */
124: SNESGetKSP(snes,&ksp);
125: KSPGetPC(ksp,&pc);
126: PCSetType(pc,PCNONE);
127: KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);
129: /*
130: Set SNES/KSP/KSP/PC runtime options, e.g.,
131: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
132: These options will override those specified above as long as
133: SNESSetFromOptions() is called _after_ any other customization
134: routines.
135: */
136: SNESSetFromOptions(snes);
138: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
139: Evaluate initial guess; then solve nonlinear system
140: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
141: if (!flg) {
142: VecSet(x,pfive);
143: } else {
144: VecGetArray(x,&xx);
145: xx[0] = 2.0; xx[1] = 3.0;
146: VecRestoreArray(x,&xx);
147: }
148: /*
149: Note: The user should initialize the vector, x, with the initial guess
150: for the nonlinear solver prior to calling SNESSolve(). In particular,
151: to employ an initial guess of zero, the user should explicitly set
152: this vector to zero by calling VecSet().
153: */
155: SNESSolve(snes,NULL,x);
156: SNESGetIterationNumber(snes,&its);
157: if (flg) {
158: Vec f;
159: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
160: SNESGetFunction(snes,&f,0,0);
161: VecView(r,PETSC_VIEWER_STDOUT_WORLD);
162: }
164: PetscPrintf(PETSC_COMM_SELF,"number of SNES iterations = %D\n\n",its);
166: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167: Free work space. All PETSc objects should be destroyed when they
168: are no longer needed.
169: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
171: VecDestroy(&x); VecDestroy(&r);
172: MatDestroy(&J); SNESDestroy(&snes);
173: return(0);
174: }
176: /* ------------------------------------------------------------------- */
179: /*
180: FormFunction1 - Evaluates nonlinear function, F(x).
182: Input Parameters:
183: . snes - the SNES context
184: . x - input vector
185: . dummy - optional user-defined context (not used here)
187: Output Parameter:
188: . f - function vector
189: */
190: static PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *dummy)
191: {
192: PetscErrorCode ierr;
193: const PetscScalar *xx;
194: PetscScalar *ff;
197: /*
198: Get pointers to vector data.
199: - For default PETSc vectors, VecGetArray() returns a pointer to
200: the data array. Otherwise, the routine is implementation dependent.
201: - You MUST call VecRestoreArray() when you no longer need access to
202: the array.
203: */
204: VecGetArrayRead(x,&xx);
205: VecGetArray(f,&ff);
207: /*
208: Compute function
209: */
210: ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0;
211: ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0;
214: /*
215: Restore vectors
216: */
217: VecRestoreArrayRead(x,&xx);
218: VecRestoreArray(f,&ff);
219: return(0);
220: }
221: /* ------------------------------------------------------------------- */
224: /*
225: FormJacobian1 - Evaluates Jacobian matrix.
227: Input Parameters:
228: . snes - the SNES context
229: . x - input vector
230: . dummy - optional user-defined context (not used here)
232: Output Parameters:
233: . jac - Jacobian matrix
234: . B - optionally different preconditioning matrix
235: . flag - flag indicating matrix structure
236: */
237: static PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
238: {
239: const PetscScalar *xx;
240: PetscScalar A[4];
241: PetscErrorCode ierr;
242: PetscInt idx[2] = {0,1};
245: /*
246: Get pointer to vector data
247: */
248: VecGetArrayRead(x,&xx);
250: /*
251: Compute Jacobian entries and insert into matrix.
252: - Since this is such a small problem, we set all entries for
253: the matrix at once.
254: */
255: A[0] = 2.0*xx[0] + xx[1]; A[1] = xx[0];
256: A[2] = xx[1]; A[3] = xx[0] + 2.0*xx[1];
257: MatSetValues(jac,2,idx,2,idx,A,INSERT_VALUES);
259: /*
260: Restore vector
261: */
262: VecRestoreArrayRead(x,&xx);
264: /*
265: Assemble matrix
266: */
267: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
268: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
269: return(0);
270: }
273: /* ------------------------------------------------------------------- */
276: static PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
277: {
278: PetscErrorCode ierr;
279: const PetscScalar *xx;
280: PetscScalar *ff;
283: /*
284: Get pointers to vector data.
285: - For default PETSc vectors, VecGetArray() returns a pointer to
286: the data array. Otherwise, the routine is implementation dependent.
287: - You MUST call VecRestoreArray() when you no longer need access to
288: the array.
289: */
290: VecGetArrayRead(x,&xx);
291: VecGetArray(f,&ff);
293: /*
294: Compute function
295: */
296: ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
297: ff[1] = xx[1];
299: /*
300: Restore vectors
301: */
302: VecRestoreArrayRead(x,&xx);
303: VecRestoreArray(f,&ff);
304: return(0);
305: }
307: /* ------------------------------------------------------------------- */
310: static PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
311: {
312: const PetscScalar *xx;
313: PetscScalar A[4];
314: PetscErrorCode ierr;
315: PetscInt idx[2] = {0,1};
318: /*
319: Get pointer to vector data
320: */
321: VecGetArrayRead(x,&xx);
323: /*
324: Compute Jacobian entries and insert into matrix.
325: - Since this is such a small problem, we set all entries for
326: the matrix at once.
327: */
328: A[0] = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
329: A[2] = 0.0; A[3] = 1.0;
330: MatSetValues(jac,2,idx,2,idx,A,INSERT_VALUES);
332: /*
333: Restore vector
334: */
335: VecRestoreArrayRead(x,&xx);
337: /*
338: Assemble matrix
339: */
340: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
341: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
342: return(0);
343: }
347: static int block_system(void)
348: {
349: SNES snes; /* nonlinear solver context */
350: KSP ksp; /* linear solver context */
351: PC pc; /* preconditioner context */
352: Vec x,r; /* solution, residual vectors */
353: Mat J; /* Jacobian matrix */
355: PetscInt its;
356: PetscScalar pfive = .5;
357: PetscBool flg;
359: Mat j11, j12, j21, j22;
360: Vec x1, x2, r1, r2;
361: Vec bv;
362: Vec bx[2];
363: Mat bA[2][2];
366: PetscPrintf(PETSC_COMM_WORLD, "\n\n========================= Block system =========================\n\n");
368: SNESCreate(PETSC_COMM_WORLD,&snes);
370: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
371: Create matrix and vector data structures; set corresponding routines
372: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
374: /*
375: Create sub vectors for solution and nonlinear function
376: */
377: VecCreateSeq(PETSC_COMM_SELF,1,&x1);
378: VecDuplicate(x1,&r1);
380: VecCreateSeq(PETSC_COMM_SELF,1,&x2);
381: VecDuplicate(x2,&r2);
383: /*
384: Create the block vectors
385: */
386: bx[0] = x1;
387: bx[1] = x2;
388: VecCreateNest(PETSC_COMM_WORLD,2,NULL,bx,&x);
389: VecAssemblyBegin(x);
390: VecAssemblyEnd(x);
391: VecDestroy(&x1);
392: VecDestroy(&x2);
394: bx[0] = r1;
395: bx[1] = r2;
396: VecCreateNest(PETSC_COMM_WORLD,2,NULL,bx,&r);
397: VecDestroy(&r1);
398: VecDestroy(&r2);
399: VecAssemblyBegin(r);
400: VecAssemblyEnd(r);
402: /*
403: Create sub Jacobian matrix data structure
404: */
405: MatCreate(PETSC_COMM_WORLD, &j11);
406: MatSetSizes(j11, 1, 1, 1, 1);
407: MatSetType(j11, MATSEQAIJ);
408: MatSetUp(j11);
410: MatCreate(PETSC_COMM_WORLD, &j12);
411: MatSetSizes(j12, 1, 1, 1, 1);
412: MatSetType(j12, MATSEQAIJ);
413: MatSetUp(j12);
415: MatCreate(PETSC_COMM_WORLD, &j21);
416: MatSetSizes(j21, 1, 1, 1, 1);
417: MatSetType(j21, MATSEQAIJ);
418: MatSetUp(j21);
420: MatCreate(PETSC_COMM_WORLD, &j22);
421: MatSetSizes(j22, PETSC_DECIDE, PETSC_DECIDE, 1, 1);
422: MatSetType(j22, MATSEQAIJ);
423: MatSetUp(j22);
424: /*
425: Create block Jacobian matrix data structure
426: */
427: bA[0][0] = j11;
428: bA[0][1] = j12;
429: bA[1][0] = j21;
430: bA[1][1] = j22;
432: MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,&bA[0][0],&J);
433: MatSetUp(J);
434: MatNestSetVecType(J,VECNEST);
435: MatDestroy(&j11);
436: MatDestroy(&j12);
437: MatDestroy(&j21);
438: MatDestroy(&j22);
440: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
441: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
443: PetscOptionsHasName(NULL,"-hard",&flg);
444: if (!flg) {
445: /*
446: Set function evaluation routine and vector.
447: */
448: SNESSetFunction(snes,r,FormFunction1_block,NULL);
450: /*
451: Set Jacobian matrix data structure and Jacobian evaluation routine
452: */
453: SNESSetJacobian(snes,J,J,FormJacobian1_block,NULL);
454: } else {
455: SNESSetFunction(snes,r,FormFunction2_block,NULL);
456: SNESSetJacobian(snes,J,J,FormJacobian2_block,NULL);
457: }
459: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
460: Customize nonlinear solver; set runtime options
461: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
463: /*
464: Set linear solver defaults for this problem. By extracting the
465: KSP, KSP, and PC contexts from the SNES context, we can then
466: directly call any KSP, KSP, and PC routines to set various options.
467: */
468: SNESGetKSP(snes,&ksp);
469: KSPGetPC(ksp,&pc);
470: PCSetType(pc,PCNONE);
471: KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);
473: /*
474: Set SNES/KSP/KSP/PC runtime options, e.g.,
475: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
476: These options will override those specified above as long as
477: SNESSetFromOptions() is called _after_ any other customization
478: routines.
479: */
480: SNESSetFromOptions(snes);
482: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
483: Evaluate initial guess; then solve nonlinear system
484: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
485: if (!flg) {
486: VecSet(x,pfive);
487: } else {
488: Vec *vecs;
489: VecNestGetSubVecs(x, NULL, &vecs);
490: bv = vecs[0];
491: /* VecBlockGetSubVec(x, 0, &bv); */
492: VecSetValue(bv, 0, 2.0, INSERT_VALUES); /* xx[0] = 2.0; */
493: VecAssemblyBegin(bv);
494: VecAssemblyEnd(bv);
496: /* VecBlockGetSubVec(x, 1, &bv); */
497: bv = vecs[1];
498: VecSetValue(bv, 0, 3.0, INSERT_VALUES); /* xx[1] = 3.0; */
499: VecAssemblyBegin(bv);
500: VecAssemblyEnd(bv);
501: }
502: /*
503: Note: The user should initialize the vector, x, with the initial guess
504: for the nonlinear solver prior to calling SNESSolve(). In particular,
505: to employ an initial guess of zero, the user should explicitly set
506: this vector to zero by calling VecSet().
507: */
508: SNESSolve(snes,NULL,x);
509: SNESGetIterationNumber(snes,&its);
510: if (flg) {
511: Vec f;
512: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
513: SNESGetFunction(snes,&f,0,0);
514: VecView(r,PETSC_VIEWER_STDOUT_WORLD);
515: }
517: PetscPrintf(PETSC_COMM_SELF,"number of SNES iterations = %D\n\n",its);
519: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
520: Free work space. All PETSc objects should be destroyed when they
521: are no longer needed.
522: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
523: VecDestroy(&x); VecDestroy(&r);
524: MatDestroy(&J); SNESDestroy(&snes);
525: return(0);
526: }
528: /* ------------------------------------------------------------------- */
531: static PetscErrorCode FormFunction1_block(SNES snes,Vec x,Vec f,void *dummy)
532: {
533: Vec *xx, *ff, x1,x2, f1,f2;
534: PetscScalar ff_0, ff_1;
535: PetscScalar xx_0, xx_1;
536: PetscInt index,nb;
540: /* get blocks for function */
541: VecNestGetSubVecs(f, &nb, &ff);
542: f1 = ff[0]; f2 = ff[1];
544: /* get blocks for solution */
545: VecNestGetSubVecs(x, &nb, &xx);
546: x1 = xx[0]; x2 = xx[1];
548: /* get solution values */
549: index = 0;
550: VecGetValues(x1,1, &index, &xx_0);
551: VecGetValues(x2,1, &index, &xx_1);
553: /* Compute function */
554: ff_0 = xx_0*xx_0 + xx_0*xx_1 - 3.0;
555: ff_1 = xx_0*xx_1 + xx_1*xx_1 - 6.0;
557: /* set function values */
558: VecSetValue(f1, index, ff_0, INSERT_VALUES);
560: VecSetValue(f2, index, ff_1, INSERT_VALUES);
562: VecAssemblyBegin(f);
563: VecAssemblyEnd(f);
564: return(0);
565: }
567: /* ------------------------------------------------------------------- */
570: static PetscErrorCode FormJacobian1_block(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
571: {
572: Vec *xx, x1,x2;
573: PetscScalar xx_0, xx_1;
574: PetscInt index,nb;
575: PetscScalar A_00, A_01, A_10, A_11;
576: Mat j11, j12, j21, j22;
577: Mat **mats;
581: /* get blocks for solution */
582: VecNestGetSubVecs(x, &nb, &xx);
583: x1 = xx[0]; x2 = xx[1];
585: /* get solution values */
586: index = 0;
587: VecGetValues(x1,1, &index, &xx_0);
588: VecGetValues(x2,1, &index, &xx_1);
590: /* get block matrices */
591: MatNestGetSubMats(jac,NULL,NULL,&mats);
592: j11 = mats[0][0];
593: j12 = mats[0][1];
594: j21 = mats[1][0];
595: j22 = mats[1][1];
597: /* compute jacobian entries */
598: A_00 = 2.0*xx_0 + xx_1;
599: A_01 = xx_0;
600: A_10 = xx_1;
601: A_11 = xx_0 + 2.0*xx_1;
603: /* set jacobian values */
604: MatSetValue(j11, 0,0, A_00, INSERT_VALUES);
605: MatSetValue(j12, 0,0, A_01, INSERT_VALUES);
606: MatSetValue(j21, 0,0, A_10, INSERT_VALUES);
607: MatSetValue(j22, 0,0, A_11, INSERT_VALUES);
609: /* Assemble sub matrix */
610: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
611: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
612: return(0);
613: }
615: /* ------------------------------------------------------------------- */
618: static PetscErrorCode FormFunction2_block(SNES snes,Vec x,Vec f,void *dummy)
619: {
620: PetscErrorCode ierr;
621: PetscScalar *ff;
622: const PetscScalar *xx;
625: /*
626: Get pointers to vector data.
627: - For default PETSc vectors, VecGetArray() returns a pointer to
628: the data array. Otherwise, the routine is implementation dependent.
629: - You MUST call VecRestoreArray() when you no longer need access to
630: the array.
631: */
632: VecGetArrayRead(x,&xx);
633: VecGetArray(f,&ff);
635: /*
636: Compute function
637: */
638: ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
639: ff[1] = xx[1];
641: /*
642: Restore vectors
643: */
644: VecRestoreArrayRead(x,&xx);
645: VecRestoreArray(f,&ff);
646: return(0);
647: }
649: /* ------------------------------------------------------------------- */
652: static PetscErrorCode FormJacobian2_block(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
653: {
654: const PetscScalar *xx;
655: PetscScalar A[4];
656: PetscErrorCode ierr;
657: PetscInt idx[2] = {0,1};
660: /*
661: Get pointer to vector data
662: */
663: VecGetArrayRead(x,&xx);
665: /*
666: Compute Jacobian entries and insert into matrix.
667: - Since this is such a small problem, we set all entries for
668: the matrix at once.
669: */
670: A[0] = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
671: A[2] = 0.0; A[3] = 1.0;
672: MatSetValues(jac,2,idx,2,idx,A,INSERT_VALUES);
674: /*
675: Restore vector
676: */
677: VecRestoreArrayRead(x,&xx);
679: /*
680: Assemble matrix
681: */
682: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
683: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
684: return(0);
685: }
689: int main(int argc,char **argv)
690: {
691: PetscMPIInt size;
694: PetscInitialize(&argc,&argv,(char*)0,help);
696: MPI_Comm_size(PETSC_COMM_WORLD,&size);
697: if (size != 1) SETERRQ(PETSC_COMM_WORLD, 1,"This is a uniprocessor example only!");
699: assembled_system();
701: block_system();
703: PetscFinalize();
704: return 0;
705: }