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: /*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/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*);
57: static PetscErrorCode assembled_system(void)
58: {
59: SNES snes; /* nonlinear solver context */
60: KSP ksp; /* linear solver context */
61: PC pc; /* preconditioner context */
62: Vec x,r; /* solution, residual vectors */
63: Mat J; /* Jacobian matrix */
64: PetscInt its;
65: PetscScalar pfive = .5,*xx;
66: PetscBool flg;
69: PetscPrintf(PETSC_COMM_WORLD, "\n\n========================= Assembled system =========================\n\n");
71: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
72: Create nonlinear solver context
73: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
75: SNESCreate(PETSC_COMM_WORLD,&snes);
77: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
78: Create matrix and vector data structures; set corresponding routines
79: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
81: /*
82: Create vectors for solution and nonlinear function
83: */
84: VecCreateSeq(PETSC_COMM_SELF,2,&x);
85: VecDuplicate(x,&r);
87: /*
88: Create Jacobian matrix data structure
89: */
90: MatCreate(PETSC_COMM_SELF,&J);
91: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);
92: MatSetFromOptions(J);
93: MatSetUp(J);
95: PetscOptionsHasName(NULL,NULL,"-hard",&flg);
96: if (!flg) {
97: /*
98: Set function evaluation routine and vector.
99: */
100: SNESSetFunction(snes,r,FormFunction1,NULL);
102: /*
103: Set Jacobian matrix data structure and Jacobian evaluation routine
104: */
105: SNESSetJacobian(snes,J,J,FormJacobian1,NULL);
106: } else {
107: SNESSetFunction(snes,r,FormFunction2,NULL);
108: SNESSetJacobian(snes,J,J,FormJacobian2,NULL);
109: }
111: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
112: Customize nonlinear solver; set runtime options
113: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
115: /*
116: Set linear solver defaults for this problem. By extracting the
117: KSP, KSP, and PC contexts from the SNES context, we can then
118: directly call any KSP, KSP, and PC routines to set various options.
119: */
120: SNESGetKSP(snes,&ksp);
121: KSPGetPC(ksp,&pc);
122: PCSetType(pc,PCNONE);
123: KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);
125: /*
126: Set SNES/KSP/KSP/PC runtime options, e.g.,
127: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
128: These options will override those specified above as long as
129: SNESSetFromOptions() is called _after_ any other customization
130: routines.
131: */
132: SNESSetFromOptions(snes);
134: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135: Evaluate initial guess; then solve nonlinear system
136: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
137: if (!flg) {
138: VecSet(x,pfive);
139: } else {
140: VecGetArray(x,&xx);
141: xx[0] = 2.0; xx[1] = 3.0;
142: VecRestoreArray(x,&xx);
143: }
144: /*
145: Note: The user should initialize the vector, x, with the initial guess
146: for the nonlinear solver prior to calling SNESSolve(). In particular,
147: to employ an initial guess of zero, the user should explicitly set
148: this vector to zero by calling VecSet().
149: */
151: SNESSolve(snes,NULL,x);
152: SNESGetIterationNumber(snes,&its);
153: if (flg) {
154: Vec f;
155: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
156: SNESGetFunction(snes,&f,0,0);
157: VecView(r,PETSC_VIEWER_STDOUT_WORLD);
158: }
159: PetscPrintf(PETSC_COMM_SELF,"number of SNES iterations = %D\n\n",its);
161: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162: Free work space. All PETSc objects should be destroyed when they
163: are no longer needed.
164: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
166: VecDestroy(&x)); PetscCall(VecDestroy(&r);
167: MatDestroy(&J)); PetscCall(SNESDestroy(&snes);
168: return 0;
169: }
171: /* ------------------------------------------------------------------- */
172: /*
173: FormFunction1 - Evaluates nonlinear function, F(x).
175: Input Parameters:
176: . snes - the SNES context
177: . x - input vector
178: . dummy - optional user-defined context (not used here)
180: Output Parameter:
181: . f - function vector
182: */
183: static PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *dummy)
184: {
185: const PetscScalar *xx;
186: PetscScalar *ff;
189: /*
190: Get pointers to vector data.
191: - For default PETSc vectors, VecGetArray() returns a pointer to
192: the data array. Otherwise, the routine is implementation dependent.
193: - You MUST call VecRestoreArray() when you no longer need access to
194: the array.
195: */
196: VecGetArrayRead(x,&xx);
197: VecGetArray(f,&ff);
199: /*
200: Compute function
201: */
202: ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0;
203: ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0;
205: /*
206: Restore vectors
207: */
208: VecRestoreArrayRead(x,&xx);
209: VecRestoreArray(f,&ff);
210: return 0;
211: }
212: /* ------------------------------------------------------------------- */
213: /*
214: FormJacobian1 - Evaluates Jacobian matrix.
216: Input Parameters:
217: . snes - the SNES context
218: . x - input vector
219: . dummy - optional user-defined context (not used here)
221: Output Parameters:
222: . jac - Jacobian matrix
223: . B - optionally different preconditioning matrix
224: . flag - flag indicating matrix structure
225: */
226: static PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
227: {
228: const PetscScalar *xx;
229: PetscScalar A[4];
230: PetscInt idx[2] = {0,1};
233: /*
234: Get pointer to vector data
235: */
236: VecGetArrayRead(x,&xx);
238: /*
239: Compute Jacobian entries and insert into matrix.
240: - Since this is such a small problem, we set all entries for
241: the matrix at once.
242: */
243: A[0] = 2.0*xx[0] + xx[1]; A[1] = xx[0];
244: A[2] = xx[1]; A[3] = xx[0] + 2.0*xx[1];
245: MatSetValues(jac,2,idx,2,idx,A,INSERT_VALUES);
247: /*
248: Restore vector
249: */
250: VecRestoreArrayRead(x,&xx);
252: /*
253: Assemble matrix
254: */
255: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
256: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
257: return 0;
258: }
260: /* ------------------------------------------------------------------- */
261: static PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
262: {
263: const PetscScalar *xx;
264: PetscScalar *ff;
267: /*
268: Get pointers to vector data.
269: - For default PETSc vectors, VecGetArray() returns a pointer to
270: the data array. Otherwise, the routine is implementation dependent.
271: - You MUST call VecRestoreArray() when you no longer need access to
272: the array.
273: */
274: VecGetArrayRead(x,&xx);
275: VecGetArray(f,&ff);
277: /*
278: Compute function
279: */
280: ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
281: ff[1] = xx[1];
283: /*
284: Restore vectors
285: */
286: VecRestoreArrayRead(x,&xx);
287: VecRestoreArray(f,&ff);
288: return 0;
289: }
291: /* ------------------------------------------------------------------- */
292: static PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
293: {
294: const PetscScalar *xx;
295: PetscScalar A[4];
296: PetscInt idx[2] = {0,1};
299: /*
300: Get pointer to vector data
301: */
302: VecGetArrayRead(x,&xx);
304: /*
305: Compute Jacobian entries and insert into matrix.
306: - Since this is such a small problem, we set all entries for
307: the matrix at once.
308: */
309: A[0] = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
310: A[2] = 0.0; A[3] = 1.0;
311: MatSetValues(jac,2,idx,2,idx,A,INSERT_VALUES);
313: /*
314: Restore vector
315: */
316: VecRestoreArrayRead(x,&xx);
318: /*
319: Assemble matrix
320: */
321: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
322: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
323: return 0;
324: }
326: static PetscErrorCode block_system(void)
327: {
328: SNES snes; /* nonlinear solver context */
329: KSP ksp; /* linear solver context */
330: PC pc; /* preconditioner context */
331: Vec x,r; /* solution, residual vectors */
332: Mat J; /* Jacobian matrix */
333: PetscInt its;
334: PetscScalar pfive = .5;
335: PetscBool flg;
337: Mat j11, j12, j21, j22;
338: Vec x1, x2, r1, r2;
339: Vec bv;
340: Vec bx[2];
341: Mat bA[2][2];
344: PetscPrintf(PETSC_COMM_WORLD, "\n\n========================= Block system =========================\n\n");
346: SNESCreate(PETSC_COMM_WORLD,&snes);
348: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
349: Create matrix and vector data structures; set corresponding routines
350: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
352: /*
353: Create sub vectors for solution and nonlinear function
354: */
355: VecCreateSeq(PETSC_COMM_SELF,1,&x1);
356: VecDuplicate(x1,&r1);
358: VecCreateSeq(PETSC_COMM_SELF,1,&x2);
359: VecDuplicate(x2,&r2);
361: /*
362: Create the block vectors
363: */
364: bx[0] = x1;
365: bx[1] = x2;
366: VecCreateNest(PETSC_COMM_WORLD,2,NULL,bx,&x);
367: VecAssemblyBegin(x);
368: VecAssemblyEnd(x);
369: VecDestroy(&x1);
370: VecDestroy(&x2);
372: bx[0] = r1;
373: bx[1] = r2;
374: VecCreateNest(PETSC_COMM_WORLD,2,NULL,bx,&r);
375: VecDestroy(&r1);
376: VecDestroy(&r2);
377: VecAssemblyBegin(r);
378: VecAssemblyEnd(r);
380: /*
381: Create sub Jacobian matrix data structure
382: */
383: MatCreate(PETSC_COMM_WORLD, &j11);
384: MatSetSizes(j11, 1, 1, 1, 1);
385: MatSetType(j11, MATSEQAIJ);
386: MatSetUp(j11);
388: MatCreate(PETSC_COMM_WORLD, &j12);
389: MatSetSizes(j12, 1, 1, 1, 1);
390: MatSetType(j12, MATSEQAIJ);
391: MatSetUp(j12);
393: MatCreate(PETSC_COMM_WORLD, &j21);
394: MatSetSizes(j21, 1, 1, 1, 1);
395: MatSetType(j21, MATSEQAIJ);
396: MatSetUp(j21);
398: MatCreate(PETSC_COMM_WORLD, &j22);
399: MatSetSizes(j22, PETSC_DECIDE, PETSC_DECIDE, 1, 1);
400: MatSetType(j22, MATSEQAIJ);
401: MatSetUp(j22);
402: /*
403: Create block Jacobian matrix data structure
404: */
405: bA[0][0] = j11;
406: bA[0][1] = j12;
407: bA[1][0] = j21;
408: bA[1][1] = j22;
410: MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,&bA[0][0],&J);
411: MatSetUp(J);
412: MatNestSetVecType(J,VECNEST);
413: MatDestroy(&j11);
414: MatDestroy(&j12);
415: MatDestroy(&j21);
416: MatDestroy(&j22);
418: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
419: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
421: PetscOptionsHasName(NULL,NULL,"-hard",&flg);
422: if (!flg) {
423: /*
424: Set function evaluation routine and vector.
425: */
426: SNESSetFunction(snes,r,FormFunction1_block,NULL);
428: /*
429: Set Jacobian matrix data structure and Jacobian evaluation routine
430: */
431: SNESSetJacobian(snes,J,J,FormJacobian1_block,NULL);
432: } else {
433: SNESSetFunction(snes,r,FormFunction2_block,NULL);
434: SNESSetJacobian(snes,J,J,FormJacobian2_block,NULL);
435: }
437: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
438: Customize nonlinear solver; set runtime options
439: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
441: /*
442: Set linear solver defaults for this problem. By extracting the
443: KSP, KSP, and PC contexts from the SNES context, we can then
444: directly call any KSP, KSP, and PC routines to set various options.
445: */
446: SNESGetKSP(snes,&ksp);
447: KSPGetPC(ksp,&pc);
448: PCSetType(pc,PCNONE);
449: KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);
451: /*
452: Set SNES/KSP/KSP/PC runtime options, e.g.,
453: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
454: These options will override those specified above as long as
455: SNESSetFromOptions() is called _after_ any other customization
456: routines.
457: */
458: SNESSetFromOptions(snes);
460: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
461: Evaluate initial guess; then solve nonlinear system
462: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
463: if (!flg) {
464: VecSet(x,pfive);
465: } else {
466: Vec *vecs;
467: VecNestGetSubVecs(x, NULL, &vecs);
468: bv = vecs[0];
469: /* VecBlockGetSubVec(x, 0, &bv); */
470: VecSetValue(bv, 0, 2.0, INSERT_VALUES); /* xx[0] = 2.0; */
471: VecAssemblyBegin(bv);
472: VecAssemblyEnd(bv);
474: /* VecBlockGetSubVec(x, 1, &bv); */
475: bv = vecs[1];
476: VecSetValue(bv, 0, 3.0, INSERT_VALUES); /* xx[1] = 3.0; */
477: VecAssemblyBegin(bv);
478: VecAssemblyEnd(bv);
479: }
480: /*
481: Note: The user should initialize the vector, x, with the initial guess
482: for the nonlinear solver prior to calling SNESSolve(). In particular,
483: to employ an initial guess of zero, the user should explicitly set
484: this vector to zero by calling VecSet().
485: */
486: SNESSolve(snes,NULL,x);
487: SNESGetIterationNumber(snes,&its);
488: if (flg) {
489: Vec f;
490: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
491: SNESGetFunction(snes,&f,0,0);
492: VecView(r,PETSC_VIEWER_STDOUT_WORLD);
493: }
495: PetscPrintf(PETSC_COMM_SELF,"number of SNES iterations = %D\n\n",its);
497: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
498: Free work space. All PETSc objects should be destroyed when they
499: are no longer needed.
500: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
501: VecDestroy(&x)); PetscCall(VecDestroy(&r);
502: MatDestroy(&J)); PetscCall(SNESDestroy(&snes);
503: return 0;
504: }
506: /* ------------------------------------------------------------------- */
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]; f2 = ff[1];
519: /* get blocks for solution */
520: VecNestGetSubVecs(x, &nb, &xx);
521: x1 = xx[0]; x2 = xx[1];
523: /* get solution values */
524: index = 0;
525: VecGetValues(x1,1, &index, &xx_0);
526: VecGetValues(x2,1, &index, &xx_1);
528: /* Compute function */
529: ff_0 = xx_0*xx_0 + xx_0*xx_1 - 3.0;
530: ff_1 = xx_0*xx_1 + xx_1*xx_1 - 6.0;
532: /* set function values */
533: VecSetValue(f1, index, ff_0, INSERT_VALUES);
535: VecSetValue(f2, index, ff_1, INSERT_VALUES);
537: VecAssemblyBegin(f);
538: VecAssemblyEnd(f);
539: return 0;
540: }
542: /* ------------------------------------------------------------------- */
543: static PetscErrorCode FormJacobian1_block(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
544: {
545: Vec *xx, x1,x2;
546: PetscScalar xx_0, xx_1;
547: PetscInt index,nb;
548: PetscScalar A_00, A_01, A_10, A_11;
549: Mat j11, j12, j21, j22;
550: Mat **mats;
553: /* get blocks for solution */
554: VecNestGetSubVecs(x, &nb, &xx);
555: x1 = xx[0]; x2 = xx[1];
557: /* get solution values */
558: index = 0;
559: VecGetValues(x1,1, &index, &xx_0);
560: VecGetValues(x2,1, &index, &xx_1);
562: /* get block matrices */
563: MatNestGetSubMats(jac,NULL,NULL,&mats);
564: j11 = mats[0][0];
565: j12 = mats[0][1];
566: j21 = mats[1][0];
567: j22 = mats[1][1];
569: /* compute jacobian entries */
570: A_00 = 2.0*xx_0 + xx_1;
571: A_01 = xx_0;
572: A_10 = xx_1;
573: A_11 = xx_0 + 2.0*xx_1;
575: /* set jacobian values */
576: MatSetValue(j11, 0,0, A_00, INSERT_VALUES);
577: MatSetValue(j12, 0,0, A_01, INSERT_VALUES);
578: MatSetValue(j21, 0,0, A_10, INSERT_VALUES);
579: MatSetValue(j22, 0,0, A_11, INSERT_VALUES);
581: /* Assemble sub matrix */
582: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
583: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
584: return 0;
585: }
587: /* ------------------------------------------------------------------- */
588: static PetscErrorCode FormFunction2_block(SNES snes,Vec x,Vec f,void *dummy)
589: {
590: PetscScalar *ff;
591: const PetscScalar *xx;
594: /*
595: Get pointers to vector data.
596: - For default PETSc vectors, VecGetArray() returns a pointer to
597: the data array. Otherwise, the routine is implementation dependent.
598: - You MUST call VecRestoreArray() when you no longer need access to
599: the array.
600: */
601: VecGetArrayRead(x,&xx);
602: VecGetArray(f,&ff);
604: /*
605: Compute function
606: */
607: ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
608: ff[1] = xx[1];
610: /*
611: Restore vectors
612: */
613: VecRestoreArrayRead(x,&xx);
614: VecRestoreArray(f,&ff);
615: return 0;
616: }
618: /* ------------------------------------------------------------------- */
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; A[1] = 0.0;
637: A[2] = 0.0; A[3] = 1.0;
638: MatSetValues(jac,2,idx,2,idx,A,INSERT_VALUES);
640: /*
641: Restore vector
642: */
643: VecRestoreArrayRead(x,&xx);
645: /*
646: Assemble matrix
647: */
648: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
649: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
650: return 0;
651: }
653: int main(int argc,char **argv)
654: {
655: PetscMPIInt size;
657: PetscInitialize(&argc,&argv,(char*)0,help);
658: MPI_Comm_size(PETSC_COMM_WORLD,&size);
660: assembled_system();
661: block_system();
662: PetscFinalize();
663: return 0;
664: }
666: /*TEST
668: test:
669: args: -snes_monitor_short
670: requires: !single
672: TEST*/