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; xx[1] = 3.0;
137: VecRestoreArray(x,&xx);
138: }
139: /*
140: Note: The user should initialize the vector, x, with the initial guess
141: for the nonlinear solver prior to calling SNESSolve(). In particular,
142: to employ an initial guess of zero, the user should explicitly set
143: this vector to zero by calling VecSet().
144: */
146: SNESSolve(snes,NULL,x);
147: SNESGetIterationNumber(snes,&its);
148: if (flg) {
149: Vec f;
150: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
151: SNESGetFunction(snes,&f,0,0);
152: VecView(r,PETSC_VIEWER_STDOUT_WORLD);
153: }
154: PetscPrintf(PETSC_COMM_SELF,"number of SNES iterations = %D\n\n",its);
156: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157: Free work space. All PETSc objects should be destroyed when they
158: are no longer needed.
159: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
161: VecDestroy(&x)); PetscCall(VecDestroy(&r);
162: MatDestroy(&J)); PetscCall(SNESDestroy(&snes);
163: return 0;
164: }
166: /* ------------------------------------------------------------------- */
167: /*
168: FormFunction1 - Evaluates nonlinear function, F(x).
170: Input Parameters:
171: . snes - the SNES context
172: . x - input vector
173: . dummy - optional user-defined context (not used here)
175: Output Parameter:
176: . f - function vector
177: */
178: static PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *dummy)
179: {
180: const PetscScalar *xx;
181: PetscScalar *ff;
184: /*
185: Get pointers to vector data.
186: - For default PETSc vectors, VecGetArray() returns a pointer to
187: the data array. Otherwise, the routine is implementation dependent.
188: - You MUST call VecRestoreArray() when you no longer need access to
189: the array.
190: */
191: VecGetArrayRead(x,&xx);
192: VecGetArray(f,&ff);
194: /*
195: Compute function
196: */
197: ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0;
198: ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0;
200: /*
201: Restore vectors
202: */
203: VecRestoreArrayRead(x,&xx);
204: VecRestoreArray(f,&ff);
205: return 0;
206: }
207: /* ------------------------------------------------------------------- */
208: /*
209: FormJacobian1 - Evaluates Jacobian matrix.
211: Input Parameters:
212: . snes - the SNES context
213: . x - input vector
214: . dummy - optional user-defined context (not used here)
216: Output Parameters:
217: . jac - Jacobian matrix
218: . B - optionally different preconditioning matrix
219: . flag - flag indicating matrix structure
220: */
221: static PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
222: {
223: const PetscScalar *xx;
224: PetscScalar A[4];
225: PetscInt idx[2] = {0,1};
228: /*
229: Get pointer to vector data
230: */
231: VecGetArrayRead(x,&xx);
233: /*
234: Compute Jacobian entries and insert into matrix.
235: - Since this is such a small problem, we set all entries for
236: the matrix at once.
237: */
238: A[0] = 2.0*xx[0] + xx[1]; A[1] = xx[0];
239: A[2] = xx[1]; A[3] = xx[0] + 2.0*xx[1];
240: MatSetValues(jac,2,idx,2,idx,A,INSERT_VALUES);
242: /*
243: Restore vector
244: */
245: VecRestoreArrayRead(x,&xx);
247: /*
248: Assemble matrix
249: */
250: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
251: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
252: return 0;
253: }
255: /* ------------------------------------------------------------------- */
256: static PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
257: {
258: const PetscScalar *xx;
259: PetscScalar *ff;
262: /*
263: Get pointers to vector data.
264: - For default PETSc vectors, VecGetArray() returns a pointer to
265: the data array. Otherwise, the routine is implementation dependent.
266: - You MUST call VecRestoreArray() when you no longer need access to
267: the array.
268: */
269: VecGetArrayRead(x,&xx);
270: VecGetArray(f,&ff);
272: /*
273: Compute function
274: */
275: ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
276: ff[1] = xx[1];
278: /*
279: Restore vectors
280: */
281: VecRestoreArrayRead(x,&xx);
282: VecRestoreArray(f,&ff);
283: return 0;
284: }
286: /* ------------------------------------------------------------------- */
287: static PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
288: {
289: const PetscScalar *xx;
290: PetscScalar A[4];
291: PetscInt idx[2] = {0,1};
294: /*
295: Get pointer to vector data
296: */
297: VecGetArrayRead(x,&xx);
299: /*
300: Compute Jacobian entries and insert into matrix.
301: - Since this is such a small problem, we set all entries for
302: the matrix at once.
303: */
304: A[0] = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
305: A[2] = 0.0; A[3] = 1.0;
306: MatSetValues(jac,2,idx,2,idx,A,INSERT_VALUES);
308: /*
309: Restore vector
310: */
311: VecRestoreArrayRead(x,&xx);
313: /*
314: Assemble matrix
315: */
316: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
317: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
318: return 0;
319: }
321: static PetscErrorCode block_system(void)
322: {
323: SNES snes; /* nonlinear solver context */
324: KSP ksp; /* linear solver context */
325: PC pc; /* preconditioner context */
326: Vec x,r; /* solution, residual vectors */
327: Mat J; /* Jacobian matrix */
328: PetscInt its;
329: PetscScalar pfive = .5;
330: PetscBool flg;
332: Mat j11, j12, j21, j22;
333: Vec x1, x2, r1, r2;
334: Vec bv;
335: Vec bx[2];
336: Mat bA[2][2];
339: PetscPrintf(PETSC_COMM_WORLD, "\n\n========================= Block system =========================\n\n");
341: SNESCreate(PETSC_COMM_WORLD,&snes);
343: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
344: Create matrix and vector data structures; set corresponding routines
345: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
347: /*
348: Create sub vectors for solution and nonlinear function
349: */
350: VecCreateSeq(PETSC_COMM_SELF,1,&x1);
351: VecDuplicate(x1,&r1);
353: VecCreateSeq(PETSC_COMM_SELF,1,&x2);
354: VecDuplicate(x2,&r2);
356: /*
357: Create the block vectors
358: */
359: bx[0] = x1;
360: bx[1] = x2;
361: VecCreateNest(PETSC_COMM_WORLD,2,NULL,bx,&x);
362: VecAssemblyBegin(x);
363: VecAssemblyEnd(x);
364: VecDestroy(&x1);
365: VecDestroy(&x2);
367: bx[0] = r1;
368: bx[1] = r2;
369: VecCreateNest(PETSC_COMM_WORLD,2,NULL,bx,&r);
370: VecDestroy(&r1);
371: VecDestroy(&r2);
372: VecAssemblyBegin(r);
373: VecAssemblyEnd(r);
375: /*
376: Create sub Jacobian matrix data structure
377: */
378: MatCreate(PETSC_COMM_WORLD, &j11);
379: MatSetSizes(j11, 1, 1, 1, 1);
380: MatSetType(j11, MATSEQAIJ);
381: MatSetUp(j11);
383: MatCreate(PETSC_COMM_WORLD, &j12);
384: MatSetSizes(j12, 1, 1, 1, 1);
385: MatSetType(j12, MATSEQAIJ);
386: MatSetUp(j12);
388: MatCreate(PETSC_COMM_WORLD, &j21);
389: MatSetSizes(j21, 1, 1, 1, 1);
390: MatSetType(j21, MATSEQAIJ);
391: MatSetUp(j21);
393: MatCreate(PETSC_COMM_WORLD, &j22);
394: MatSetSizes(j22, PETSC_DECIDE, PETSC_DECIDE, 1, 1);
395: MatSetType(j22, MATSEQAIJ);
396: MatSetUp(j22);
397: /*
398: Create block Jacobian matrix data structure
399: */
400: bA[0][0] = j11;
401: bA[0][1] = j12;
402: bA[1][0] = j21;
403: bA[1][1] = j22;
405: MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,&bA[0][0],&J);
406: MatSetUp(J);
407: MatNestSetVecType(J,VECNEST);
408: MatDestroy(&j11);
409: MatDestroy(&j12);
410: MatDestroy(&j21);
411: MatDestroy(&j22);
413: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
414: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
416: PetscOptionsHasName(NULL,NULL,"-hard",&flg);
417: if (!flg) {
418: /*
419: Set function evaluation routine and vector.
420: */
421: SNESSetFunction(snes,r,FormFunction1_block,NULL);
423: /*
424: Set Jacobian matrix data structure and Jacobian evaluation routine
425: */
426: SNESSetJacobian(snes,J,J,FormJacobian1_block,NULL);
427: } else {
428: SNESSetFunction(snes,r,FormFunction2_block,NULL);
429: SNESSetJacobian(snes,J,J,FormJacobian2_block,NULL);
430: }
432: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
433: Customize nonlinear solver; set runtime options
434: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
436: /*
437: Set linear solver defaults for this problem. By extracting the
438: KSP, KSP, and PC contexts from the SNES context, we can then
439: directly call any KSP, KSP, and PC routines to set various options.
440: */
441: SNESGetKSP(snes,&ksp);
442: KSPGetPC(ksp,&pc);
443: PCSetType(pc,PCNONE);
444: KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);
446: /*
447: Set SNES/KSP/KSP/PC runtime options, e.g.,
448: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
449: These options will override those specified above as long as
450: SNESSetFromOptions() is called _after_ any other customization
451: routines.
452: */
453: SNESSetFromOptions(snes);
455: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
456: Evaluate initial guess; then solve nonlinear system
457: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
458: if (!flg) {
459: VecSet(x,pfive);
460: } else {
461: Vec *vecs;
462: VecNestGetSubVecs(x, NULL, &vecs);
463: bv = vecs[0];
464: /* VecBlockGetSubVec(x, 0, &bv); */
465: VecSetValue(bv, 0, 2.0, INSERT_VALUES); /* xx[0] = 2.0; */
466: VecAssemblyBegin(bv);
467: VecAssemblyEnd(bv);
469: /* VecBlockGetSubVec(x, 1, &bv); */
470: bv = vecs[1];
471: VecSetValue(bv, 0, 3.0, INSERT_VALUES); /* xx[1] = 3.0; */
472: VecAssemblyBegin(bv);
473: VecAssemblyEnd(bv);
474: }
475: /*
476: Note: The user should initialize the vector, x, with the initial guess
477: for the nonlinear solver prior to calling SNESSolve(). In particular,
478: to employ an initial guess of zero, the user should explicitly set
479: this vector to zero by calling VecSet().
480: */
481: SNESSolve(snes,NULL,x);
482: SNESGetIterationNumber(snes,&its);
483: if (flg) {
484: Vec f;
485: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
486: SNESGetFunction(snes,&f,0,0);
487: VecView(r,PETSC_VIEWER_STDOUT_WORLD);
488: }
490: PetscPrintf(PETSC_COMM_SELF,"number of SNES iterations = %D\n\n",its);
492: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
493: Free work space. All PETSc objects should be destroyed when they
494: are no longer needed.
495: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
496: VecDestroy(&x)); PetscCall(VecDestroy(&r);
497: MatDestroy(&J)); PetscCall(SNESDestroy(&snes);
498: return 0;
499: }
501: /* ------------------------------------------------------------------- */
502: static PetscErrorCode FormFunction1_block(SNES snes,Vec x,Vec f,void *dummy)
503: {
504: Vec *xx, *ff, x1,x2, f1,f2;
505: PetscScalar ff_0, ff_1;
506: PetscScalar xx_0, xx_1;
507: PetscInt index,nb;
510: /* get blocks for function */
511: VecNestGetSubVecs(f, &nb, &ff);
512: f1 = ff[0]; f2 = ff[1];
514: /* get blocks for solution */
515: VecNestGetSubVecs(x, &nb, &xx);
516: x1 = xx[0]; x2 = xx[1];
518: /* get solution values */
519: index = 0;
520: VecGetValues(x1,1, &index, &xx_0);
521: VecGetValues(x2,1, &index, &xx_1);
523: /* Compute function */
524: ff_0 = xx_0*xx_0 + xx_0*xx_1 - 3.0;
525: ff_1 = xx_0*xx_1 + xx_1*xx_1 - 6.0;
527: /* set function values */
528: VecSetValue(f1, index, ff_0, INSERT_VALUES);
530: VecSetValue(f2, index, ff_1, INSERT_VALUES);
532: VecAssemblyBegin(f);
533: VecAssemblyEnd(f);
534: return 0;
535: }
537: /* ------------------------------------------------------------------- */
538: static PetscErrorCode FormJacobian1_block(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
539: {
540: Vec *xx, x1,x2;
541: PetscScalar xx_0, xx_1;
542: PetscInt index,nb;
543: PetscScalar A_00, A_01, A_10, A_11;
544: Mat j11, j12, j21, j22;
545: Mat **mats;
548: /* get blocks for solution */
549: VecNestGetSubVecs(x, &nb, &xx);
550: x1 = xx[0]; x2 = xx[1];
552: /* get solution values */
553: index = 0;
554: VecGetValues(x1,1, &index, &xx_0);
555: VecGetValues(x2,1, &index, &xx_1);
557: /* get block matrices */
558: MatNestGetSubMats(jac,NULL,NULL,&mats);
559: j11 = mats[0][0];
560: j12 = mats[0][1];
561: j21 = mats[1][0];
562: j22 = mats[1][1];
564: /* compute jacobian entries */
565: A_00 = 2.0*xx_0 + xx_1;
566: A_01 = xx_0;
567: A_10 = xx_1;
568: A_11 = xx_0 + 2.0*xx_1;
570: /* set jacobian values */
571: MatSetValue(j11, 0,0, A_00, INSERT_VALUES);
572: MatSetValue(j12, 0,0, A_01, INSERT_VALUES);
573: MatSetValue(j21, 0,0, A_10, INSERT_VALUES);
574: MatSetValue(j22, 0,0, A_11, INSERT_VALUES);
576: /* Assemble sub matrix */
577: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
578: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
579: return 0;
580: }
582: /* ------------------------------------------------------------------- */
583: static PetscErrorCode FormFunction2_block(SNES snes,Vec x,Vec f,void *dummy)
584: {
585: PetscScalar *ff;
586: const PetscScalar *xx;
589: /*
590: Get pointers to vector data.
591: - For default PETSc vectors, VecGetArray() returns a pointer to
592: the data array. Otherwise, the routine is implementation dependent.
593: - You MUST call VecRestoreArray() when you no longer need access to
594: the array.
595: */
596: VecGetArrayRead(x,&xx);
597: VecGetArray(f,&ff);
599: /*
600: Compute function
601: */
602: ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
603: ff[1] = xx[1];
605: /*
606: Restore vectors
607: */
608: VecRestoreArrayRead(x,&xx);
609: VecRestoreArray(f,&ff);
610: return 0;
611: }
613: /* ------------------------------------------------------------------- */
614: static PetscErrorCode FormJacobian2_block(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
615: {
616: const PetscScalar *xx;
617: PetscScalar A[4];
618: PetscInt idx[2] = {0,1};
621: /*
622: Get pointer to vector data
623: */
624: VecGetArrayRead(x,&xx);
626: /*
627: Compute Jacobian entries and insert into matrix.
628: - Since this is such a small problem, we set all entries for
629: the matrix at once.
630: */
631: A[0] = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
632: A[2] = 0.0; A[3] = 1.0;
633: MatSetValues(jac,2,idx,2,idx,A,INSERT_VALUES);
635: /*
636: Restore vector
637: */
638: VecRestoreArrayRead(x,&xx);
640: /*
641: Assemble matrix
642: */
643: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
644: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
645: return 0;
646: }
648: int main(int argc,char **argv)
649: {
650: PetscMPIInt size;
652: PetscInitialize(&argc,&argv,(char*)0,help);
653: MPI_Comm_size(PETSC_COMM_WORLD,&size);
655: assembled_system();
656: block_system();
657: PetscFinalize();
658: return 0;
659: }
661: /*TEST
663: test:
664: args: -snes_monitor_short
665: requires: !single
667: TEST*/