Actual source code: ex17.c
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: /*
5: Include "petscsnes.h" so that we can use SNES solvers. Note that this
6: file automatically includes:
7: petscsys.h - base PETSc routines petscvec.h - vectors
8: petscsys.h - system routines petscmat.h - matrices
9: petscis.h - index sets petscksp.h - Krylov subspace methods
10: petscviewer.h - viewers petscpc.h - preconditioners
11: petscksp.h - linear solvers
12: */
13: #include <petscsnes.h>
15: /*
16: This example is block version of the test found at
17: ${PETSC_DIR}/src/snes/tutorials/ex1.c
18: In this test we replace the Jacobian systems
19: [J]{x} = {F}
20: where
22: [J] = (j_00, j_01), {x} = (x_0, x_1)^T, {F} = (f_0, f_1)^T
23: (j_10, j_11)
24: where [J] \in \mathbb^{2 \times 2}, {x},{F} \in \mathbb^{2 \times 1},
26: with a block system in which each block is of length 1.
27: i.e. The block system is thus
29: [J] = ([j00], [j01]), {x} = ({x0}, {x1})^T, {F} = ({f0}, {f1})^T
30: ([j10], [j11])
31: where
32: [j00], [j01], [j10], [j11] \in \mathbb^{1 \times 1}
33: {x0}, {x1}, {f0}, {f1} \in \mathbb^{1 \times 1}
35: In practice we would not bother defing blocks of size one, and would instead assemble the
36: full system. This is just a simple test to illustrate how to manipulate the blocks and
37: to confirm the implementation is correct.
38: */
40: /*
41: User-defined routines
42: */
43: static PetscErrorCode FormJacobian1(SNES, Vec, Mat, Mat, void *);
44: static PetscErrorCode FormFunction1(SNES, Vec, Vec, void *);
45: static PetscErrorCode FormJacobian2(SNES, Vec, Mat, Mat, void *);
46: static PetscErrorCode FormFunction2(SNES, Vec, Vec, void *);
47: static PetscErrorCode FormJacobian1_block(SNES, Vec, Mat, Mat, void *);
48: static PetscErrorCode FormFunction1_block(SNES, Vec, Vec, void *);
49: static PetscErrorCode FormJacobian2_block(SNES, Vec, Mat, Mat, void *);
50: static PetscErrorCode FormFunction2_block(SNES, Vec, Vec, void *);
52: static PetscErrorCode assembled_system(void)
53: {
54: SNES snes; /* nonlinear solver context */
55: KSP ksp; /* linear solver context */
56: PC pc; /* preconditioner context */
57: Vec x, r; /* solution, residual vectors */
58: Mat J; /* Jacobian matrix */
59: PetscInt its;
60: PetscScalar pfive = .5, *xx;
61: PetscBool flg;
64: PetscPrintf(PETSC_COMM_WORLD, "\n\n========================= Assembled system =========================\n\n");
66: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
67: Create nonlinear solver context
68: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
70: SNESCreate(PETSC_COMM_WORLD, &snes);
72: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
73: Create matrix and vector data structures; set corresponding routines
74: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
76: /*
77: Create vectors for solution and nonlinear function
78: */
79: VecCreateSeq(PETSC_COMM_SELF, 2, &x);
80: VecDuplicate(x, &r);
82: /*
83: Create Jacobian matrix data structure
84: */
85: MatCreate(PETSC_COMM_SELF, &J);
86: MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, 2, 2);
87: MatSetFromOptions(J);
88: MatSetUp(J);
90: PetscOptionsHasName(NULL, NULL, "-hard", &flg);
91: if (!flg) {
92: /*
93: Set function evaluation routine and vector.
94: */
95: SNESSetFunction(snes, r, FormFunction1, NULL);
97: /*
98: Set Jacobian matrix data structure and Jacobian evaluation routine
99: */
100: SNESSetJacobian(snes, J, J, FormJacobian1, NULL);
101: } else {
102: SNESSetFunction(snes, r, FormFunction2, NULL);
103: SNESSetJacobian(snes, J, J, FormJacobian2, NULL);
104: }
106: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: Customize nonlinear solver; set runtime options
108: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110: /*
111: Set linear solver defaults for this problem. By extracting the
112: KSP, KSP, and PC contexts from the SNES context, we can then
113: directly call any KSP, KSP, and PC routines to set various options.
114: */
115: SNESGetKSP(snes, &ksp);
116: KSPGetPC(ksp, &pc);
117: PCSetType(pc, PCNONE);
118: KSPSetTolerances(ksp, 1.e-4, PETSC_DEFAULT, PETSC_DEFAULT, 20);
120: /*
121: Set SNES/KSP/KSP/PC runtime options, e.g.,
122: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
123: These options will override those specified above as long as
124: SNESSetFromOptions() is called _after_ any other customization
125: routines.
126: */
127: SNESSetFromOptions(snes);
129: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130: Evaluate initial guess; then solve nonlinear system
131: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
132: if (!flg) {
133: VecSet(x, pfive);
134: } else {
135: VecGetArray(x, &xx);
136: xx[0] = 2.0;
137: xx[1] = 3.0;
138: VecRestoreArray(x, &xx);
139: }
140: /*
141: Note: The user should initialize the vector, x, with the initial guess
142: for the nonlinear solver prior to calling SNESSolve(). In particular,
143: to employ an initial guess of zero, the user should explicitly set
144: this vector to zero by calling VecSet().
145: */
147: SNESSolve(snes, NULL, x);
148: SNESGetIterationNumber(snes, &its);
149: if (flg) {
150: Vec f;
151: VecView(x, PETSC_VIEWER_STDOUT_WORLD);
152: SNESGetFunction(snes, &f, 0, 0);
153: VecView(r, PETSC_VIEWER_STDOUT_WORLD);
154: }
155: PetscPrintf(PETSC_COMM_SELF, "number of SNES iterations = %" PetscInt_FMT "\n\n", its);
157: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158: Free work space. All PETSc objects should be destroyed when they
159: are no longer needed.
160: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
162: VecDestroy(&x);
163: VecDestroy(&r);
164: MatDestroy(&J);
165: SNESDestroy(&snes);
166: return 0;
167: }
169: /*
170: FormFunction1 - Evaluates nonlinear function, F(x).
172: Input Parameters:
173: . snes - the SNES context
174: . x - input vector
175: . dummy - optional user-defined context (not used here)
177: Output Parameter:
178: . f - function vector
179: */
180: static PetscErrorCode FormFunction1(SNES snes, Vec x, Vec f, void *dummy)
181: {
182: const PetscScalar *xx;
183: PetscScalar *ff;
186: /*
187: Get pointers to vector data.
188: - For default PETSc vectors, VecGetArray() returns a pointer to
189: the data array. Otherwise, the routine is implementation dependent.
190: - You MUST call VecRestoreArray() when you no longer need access to
191: the array.
192: */
193: VecGetArrayRead(x, &xx);
194: VecGetArray(f, &ff);
196: /*
197: Compute function
198: */
199: ff[0] = xx[0] * xx[0] + xx[0] * xx[1] - 3.0;
200: ff[1] = xx[0] * xx[1] + xx[1] * xx[1] - 6.0;
202: /*
203: Restore vectors
204: */
205: VecRestoreArrayRead(x, &xx);
206: VecRestoreArray(f, &ff);
207: return 0;
208: }
210: /*
211: FormJacobian1 - Evaluates Jacobian matrix.
213: Input Parameters:
214: . snes - the SNES context
215: . x - input vector
216: . dummy - optional user-defined context (not used here)
218: Output Parameters:
219: . jac - Jacobian matrix
220: . B - optionally different preconditioning matrix
221: . flag - flag indicating matrix structure
222: */
223: static PetscErrorCode FormJacobian1(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
224: {
225: const PetscScalar *xx;
226: PetscScalar A[4];
227: PetscInt idx[2] = {0, 1};
230: /*
231: Get pointer to vector data
232: */
233: VecGetArrayRead(x, &xx);
235: /*
236: Compute Jacobian entries and insert into matrix.
237: - Since this is such a small problem, we set all entries for
238: the matrix at once.
239: */
240: A[0] = 2.0 * xx[0] + xx[1];
241: A[1] = xx[0];
242: A[2] = xx[1];
243: A[3] = xx[0] + 2.0 * xx[1];
244: MatSetValues(jac, 2, idx, 2, idx, A, INSERT_VALUES);
246: /*
247: Restore vector
248: */
249: VecRestoreArrayRead(x, &xx);
251: /*
252: Assemble matrix
253: */
254: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
255: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
256: return 0;
257: }
259: static PetscErrorCode FormFunction2(SNES snes, Vec x, Vec f, void *dummy)
260: {
261: const PetscScalar *xx;
262: PetscScalar *ff;
265: /*
266: Get pointers to vector data.
267: - For default PETSc vectors, VecGetArray() returns a pointer to
268: the data array. Otherwise, the routine is implementation dependent.
269: - You MUST call VecRestoreArray() when you no longer need access to
270: the array.
271: */
272: VecGetArrayRead(x, &xx);
273: VecGetArray(f, &ff);
275: /*
276: Compute function
277: */
278: ff[0] = PetscSinScalar(3.0 * xx[0]) + xx[0];
279: ff[1] = xx[1];
281: /*
282: Restore vectors
283: */
284: VecRestoreArrayRead(x, &xx);
285: VecRestoreArray(f, &ff);
286: return 0;
287: }
289: static PetscErrorCode FormJacobian2(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
290: {
291: const PetscScalar *xx;
292: PetscScalar A[4];
293: PetscInt idx[2] = {0, 1};
296: /*
297: Get pointer to vector data
298: */
299: VecGetArrayRead(x, &xx);
301: /*
302: Compute Jacobian entries and insert into matrix.
303: - Since this is such a small problem, we set all entries for
304: the matrix at once.
305: */
306: A[0] = 3.0 * PetscCosScalar(3.0 * xx[0]) + 1.0;
307: A[1] = 0.0;
308: A[2] = 0.0;
309: A[3] = 1.0;
310: MatSetValues(jac, 2, idx, 2, idx, A, INSERT_VALUES);
312: /*
313: Restore vector
314: */
315: VecRestoreArrayRead(x, &xx);
317: /*
318: Assemble matrix
319: */
320: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
321: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
322: return 0;
323: }
325: static PetscErrorCode block_system(void)
326: {
327: SNES snes; /* nonlinear solver context */
328: KSP ksp; /* linear solver context */
329: PC pc; /* preconditioner context */
330: Vec x, r; /* solution, residual vectors */
331: Mat J; /* Jacobian matrix */
332: PetscInt its;
333: PetscScalar pfive = .5;
334: PetscBool flg;
336: Mat j11, j12, j21, j22;
337: Vec x1, x2, r1, r2;
338: Vec bv;
339: Vec bx[2];
340: Mat bA[2][2];
343: PetscPrintf(PETSC_COMM_WORLD, "\n\n========================= Block system =========================\n\n");
345: SNESCreate(PETSC_COMM_WORLD, &snes);
347: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
348: Create matrix and vector data structures; set corresponding routines
349: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
351: /*
352: Create sub vectors for solution and nonlinear function
353: */
354: VecCreateSeq(PETSC_COMM_SELF, 1, &x1);
355: VecDuplicate(x1, &r1);
357: VecCreateSeq(PETSC_COMM_SELF, 1, &x2);
358: VecDuplicate(x2, &r2);
360: /*
361: Create the block vectors
362: */
363: bx[0] = x1;
364: bx[1] = x2;
365: VecCreateNest(PETSC_COMM_WORLD, 2, NULL, bx, &x);
366: VecAssemblyBegin(x);
367: VecAssemblyEnd(x);
368: VecDestroy(&x1);
369: VecDestroy(&x2);
371: bx[0] = r1;
372: bx[1] = r2;
373: VecCreateNest(PETSC_COMM_WORLD, 2, NULL, bx, &r);
374: VecDestroy(&r1);
375: VecDestroy(&r2);
376: VecAssemblyBegin(r);
377: VecAssemblyEnd(r);
379: /*
380: Create sub Jacobian matrix data structure
381: */
382: MatCreate(PETSC_COMM_WORLD, &j11);
383: MatSetSizes(j11, 1, 1, 1, 1);
384: MatSetType(j11, MATSEQAIJ);
385: MatSetUp(j11);
387: MatCreate(PETSC_COMM_WORLD, &j12);
388: MatSetSizes(j12, 1, 1, 1, 1);
389: MatSetType(j12, MATSEQAIJ);
390: MatSetUp(j12);
392: MatCreate(PETSC_COMM_WORLD, &j21);
393: MatSetSizes(j21, 1, 1, 1, 1);
394: MatSetType(j21, MATSEQAIJ);
395: MatSetUp(j21);
397: MatCreate(PETSC_COMM_WORLD, &j22);
398: MatSetSizes(j22, PETSC_DECIDE, PETSC_DECIDE, 1, 1);
399: MatSetType(j22, MATSEQAIJ);
400: MatSetUp(j22);
401: /*
402: Create block Jacobian matrix data structure
403: */
404: bA[0][0] = j11;
405: bA[0][1] = j12;
406: bA[1][0] = j21;
407: bA[1][1] = j22;
409: MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, &bA[0][0], &J);
410: MatSetUp(J);
411: MatNestSetVecType(J, VECNEST);
412: MatDestroy(&j11);
413: MatDestroy(&j12);
414: MatDestroy(&j21);
415: MatDestroy(&j22);
417: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
418: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
420: PetscOptionsHasName(NULL, NULL, "-hard", &flg);
421: if (!flg) {
422: /*
423: Set function evaluation routine and vector.
424: */
425: SNESSetFunction(snes, r, FormFunction1_block, NULL);
427: /*
428: Set Jacobian matrix data structure and Jacobian evaluation routine
429: */
430: SNESSetJacobian(snes, J, J, FormJacobian1_block, NULL);
431: } else {
432: SNESSetFunction(snes, r, FormFunction2_block, NULL);
433: SNESSetJacobian(snes, J, J, FormJacobian2_block, NULL);
434: }
436: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
437: Customize nonlinear solver; set runtime options
438: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
440: /*
441: Set linear solver defaults for this problem. By extracting the
442: KSP, KSP, and PC contexts from the SNES context, we can then
443: directly call any KSP, KSP, and PC routines to set various options.
444: */
445: SNESGetKSP(snes, &ksp);
446: KSPGetPC(ksp, &pc);
447: PCSetType(pc, PCNONE);
448: KSPSetTolerances(ksp, 1.e-4, PETSC_DEFAULT, PETSC_DEFAULT, 20);
450: /*
451: Set SNES/KSP/KSP/PC runtime options, e.g.,
452: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
453: These options will override those specified above as long as
454: SNESSetFromOptions() is called _after_ any other customization
455: routines.
456: */
457: SNESSetFromOptions(snes);
459: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
460: Evaluate initial guess; then solve nonlinear system
461: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
462: if (!flg) {
463: VecSet(x, pfive);
464: } else {
465: Vec *vecs;
466: VecNestGetSubVecs(x, NULL, &vecs);
467: bv = vecs[0];
468: /* VecBlockGetSubVec(x, 0, &bv); */
469: VecSetValue(bv, 0, 2.0, INSERT_VALUES); /* xx[0] = 2.0; */
470: VecAssemblyBegin(bv);
471: VecAssemblyEnd(bv);
473: /* VecBlockGetSubVec(x, 1, &bv); */
474: bv = vecs[1];
475: VecSetValue(bv, 0, 3.0, INSERT_VALUES); /* xx[1] = 3.0; */
476: VecAssemblyBegin(bv);
477: VecAssemblyEnd(bv);
478: }
479: /*
480: Note: The user should initialize the vector, x, with the initial guess
481: for the nonlinear solver prior to calling SNESSolve(). In particular,
482: to employ an initial guess of zero, the user should explicitly set
483: this vector to zero by calling VecSet().
484: */
485: SNESSolve(snes, NULL, x);
486: SNESGetIterationNumber(snes, &its);
487: if (flg) {
488: Vec f;
489: VecView(x, PETSC_VIEWER_STDOUT_WORLD);
490: SNESGetFunction(snes, &f, 0, 0);
491: VecView(r, PETSC_VIEWER_STDOUT_WORLD);
492: }
494: PetscPrintf(PETSC_COMM_SELF, "number of SNES iterations = %" PetscInt_FMT "\n\n", its);
496: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
497: Free work space. All PETSc objects should be destroyed when they
498: are no longer needed.
499: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
500: VecDestroy(&x);
501: VecDestroy(&r);
502: MatDestroy(&J);
503: SNESDestroy(&snes);
504: return 0;
505: }
507: static PetscErrorCode FormFunction1_block(SNES snes, Vec x, Vec f, void *dummy)
508: {
509: Vec *xx, *ff, x1, x2, f1, f2;
510: PetscScalar ff_0, ff_1;
511: PetscScalar xx_0, xx_1;
512: PetscInt index, nb;
515: /* get blocks for function */
516: VecNestGetSubVecs(f, &nb, &ff);
517: f1 = ff[0];
518: f2 = ff[1];
520: /* get blocks for solution */
521: VecNestGetSubVecs(x, &nb, &xx);
522: x1 = xx[0];
523: x2 = xx[1];
525: /* get solution values */
526: index = 0;
527: VecGetValues(x1, 1, &index, &xx_0);
528: VecGetValues(x2, 1, &index, &xx_1);
530: /* Compute function */
531: ff_0 = xx_0 * xx_0 + xx_0 * xx_1 - 3.0;
532: ff_1 = xx_0 * xx_1 + xx_1 * xx_1 - 6.0;
534: /* set function values */
535: VecSetValue(f1, index, ff_0, INSERT_VALUES);
537: VecSetValue(f2, index, ff_1, INSERT_VALUES);
539: VecAssemblyBegin(f);
540: VecAssemblyEnd(f);
541: return 0;
542: }
544: static PetscErrorCode FormJacobian1_block(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
545: {
546: Vec *xx, x1, x2;
547: PetscScalar xx_0, xx_1;
548: PetscInt index, nb;
549: PetscScalar A_00, A_01, A_10, A_11;
550: Mat j11, j12, j21, j22;
551: Mat **mats;
554: /* get blocks for solution */
555: VecNestGetSubVecs(x, &nb, &xx);
556: x1 = xx[0];
557: x2 = xx[1];
559: /* get solution values */
560: index = 0;
561: VecGetValues(x1, 1, &index, &xx_0);
562: VecGetValues(x2, 1, &index, &xx_1);
564: /* get block matrices */
565: MatNestGetSubMats(jac, NULL, NULL, &mats);
566: j11 = mats[0][0];
567: j12 = mats[0][1];
568: j21 = mats[1][0];
569: j22 = mats[1][1];
571: /* compute jacobian entries */
572: A_00 = 2.0 * xx_0 + xx_1;
573: A_01 = xx_0;
574: A_10 = xx_1;
575: A_11 = xx_0 + 2.0 * xx_1;
577: /* set jacobian values */
578: MatSetValue(j11, 0, 0, A_00, INSERT_VALUES);
579: MatSetValue(j12, 0, 0, A_01, INSERT_VALUES);
580: MatSetValue(j21, 0, 0, A_10, INSERT_VALUES);
581: MatSetValue(j22, 0, 0, A_11, INSERT_VALUES);
583: /* Assemble sub matrix */
584: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
585: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
586: return 0;
587: }
589: static PetscErrorCode FormFunction2_block(SNES snes, Vec x, Vec f, void *dummy)
590: {
591: PetscScalar *ff;
592: const PetscScalar *xx;
595: /*
596: Get pointers to vector data.
597: - For default PETSc vectors, VecGetArray() returns a pointer to
598: the data array. Otherwise, the routine is implementation dependent.
599: - You MUST call VecRestoreArray() when you no longer need access to
600: the array.
601: */
602: VecGetArrayRead(x, &xx);
603: VecGetArray(f, &ff);
605: /*
606: Compute function
607: */
608: ff[0] = PetscSinScalar(3.0 * xx[0]) + xx[0];
609: ff[1] = xx[1];
611: /*
612: Restore vectors
613: */
614: VecRestoreArrayRead(x, &xx);
615: VecRestoreArray(f, &ff);
616: return 0;
617: }
619: static PetscErrorCode FormJacobian2_block(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
620: {
621: const PetscScalar *xx;
622: PetscScalar A[4];
623: PetscInt idx[2] = {0, 1};
626: /*
627: Get pointer to vector data
628: */
629: VecGetArrayRead(x, &xx);
631: /*
632: Compute Jacobian entries and insert into matrix.
633: - Since this is such a small problem, we set all entries for
634: the matrix at once.
635: */
636: A[0] = 3.0 * PetscCosScalar(3.0 * xx[0]) + 1.0;
637: A[1] = 0.0;
638: A[2] = 0.0;
639: A[3] = 1.0;
640: MatSetValues(jac, 2, idx, 2, idx, A, INSERT_VALUES);
642: /*
643: Restore vector
644: */
645: VecRestoreArrayRead(x, &xx);
647: /*
648: Assemble matrix
649: */
650: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
651: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
652: return 0;
653: }
655: int main(int argc, char **argv)
656: {
657: PetscMPIInt size;
660: PetscInitialize(&argc, &argv, (char *)0, help);
661: MPI_Comm_size(PETSC_COMM_WORLD, &size);
663: assembled_system();
664: block_system();
665: PetscFinalize();
666: return 0;
667: }
669: /*TEST
671: test:
672: args: -snes_monitor_short
673: requires: !single
675: TEST*/