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";
3
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>
14
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
21
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},
25
26 with a block system in which each block is of length 1.
27 i.e. The block system is thus
28
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}
34
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 */
39
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 *);
51
assembled_system(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;
62
63 PetscFunctionBeginUser;
64 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n\n========================= Assembled system =========================\n\n"));
65
66 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
67 Create nonlinear solver context
68 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
69
70 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
71
72 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
73 Create matrix and vector data structures; set corresponding routines
74 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
75
76 /*
77 Create vectors for solution and nonlinear function
78 */
79 PetscCall(VecCreateSeq(PETSC_COMM_SELF, 2, &x));
80 PetscCall(VecDuplicate(x, &r));
81
82 /*
83 Create Jacobian matrix data structure
84 */
85 PetscCall(MatCreate(PETSC_COMM_SELF, &J));
86 PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, 2, 2));
87 PetscCall(MatSetFromOptions(J));
88 PetscCall(MatSetUp(J));
89
90 PetscCall(PetscOptionsHasName(NULL, NULL, "-hard", &flg));
91 if (!flg) {
92 /*
93 Set function evaluation routine and vector.
94 */
95 PetscCall(SNESSetFunction(snes, r, FormFunction1, NULL));
96
97 /*
98 Set Jacobian matrix data structure and Jacobian evaluation routine
99 */
100 PetscCall(SNESSetJacobian(snes, J, J, FormJacobian1, NULL));
101 } else {
102 PetscCall(SNESSetFunction(snes, r, FormFunction2, NULL));
103 PetscCall(SNESSetJacobian(snes, J, J, FormJacobian2, NULL));
104 }
105
106 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107 Customize nonlinear solver; set runtime options
108 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109
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 PetscCall(SNESGetKSP(snes, &ksp));
116 PetscCall(KSPGetPC(ksp, &pc));
117 PetscCall(PCSetType(pc, PCNONE));
118 PetscCall(KSPSetTolerances(ksp, 1.e-4, PETSC_CURRENT, PETSC_CURRENT, 20));
119
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 PetscCall(SNESSetFromOptions(snes));
128
129 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130 Evaluate initial guess; then solve nonlinear system
131 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
132 if (!flg) {
133 PetscCall(VecSet(x, pfive));
134 } else {
135 PetscCall(VecGetArray(x, &xx));
136 xx[0] = 2.0;
137 xx[1] = 3.0;
138 PetscCall(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 */
146
147 PetscCall(SNESSolve(snes, NULL, x));
148 PetscCall(SNESGetIterationNumber(snes, &its));
149 if (flg) {
150 Vec f;
151 PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
152 PetscCall(SNESGetFunction(snes, &f, 0, 0));
153 PetscCall(VecView(r, PETSC_VIEWER_STDOUT_WORLD));
154 }
155 PetscCall(PetscPrintf(PETSC_COMM_SELF, "number of SNES iterations = %" PetscInt_FMT "\n\n", its));
156
157 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158 Free work space. All PETSc objects should be destroyed when they
159 are no longer needed.
160 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
161
162 PetscCall(VecDestroy(&x));
163 PetscCall(VecDestroy(&r));
164 PetscCall(MatDestroy(&J));
165 PetscCall(SNESDestroy(&snes));
166 PetscFunctionReturn(PETSC_SUCCESS);
167 }
168
169 /*
170 FormFunction1 - Evaluates nonlinear function, F(x).
171
172 Input Parameters:
173 . snes - the SNES context
174 . x - input vector
175 . dummy - optional user-defined context (not used here)
176
177 Output Parameter:
178 . f - function vector
179 */
FormFunction1(SNES snes,Vec x,Vec f,void * dummy)180 static PetscErrorCode FormFunction1(SNES snes, Vec x, Vec f, void *dummy)
181 {
182 const PetscScalar *xx;
183 PetscScalar *ff;
184
185 PetscFunctionBeginUser;
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 PetscCall(VecGetArrayRead(x, &xx));
194 PetscCall(VecGetArray(f, &ff));
195
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;
201
202 /*
203 Restore vectors
204 */
205 PetscCall(VecRestoreArrayRead(x, &xx));
206 PetscCall(VecRestoreArray(f, &ff));
207 PetscFunctionReturn(PETSC_SUCCESS);
208 }
209
210 /*
211 FormJacobian1 - Evaluates Jacobian matrix.
212
213 Input Parameters:
214 . snes - the SNES context
215 . x - input vector
216 . dummy - optional user-defined context (not used here)
217
218 Output Parameters:
219 . jac - Jacobian matrix
220 . B - optionally different matrix used to construct the preconditioner
221
222 */
FormJacobian1(SNES snes,Vec x,Mat jac,Mat B,void * dummy)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};
228
229 PetscFunctionBeginUser;
230 /*
231 Get pointer to vector data
232 */
233 PetscCall(VecGetArrayRead(x, &xx));
234
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 PetscCall(MatSetValues(jac, 2, idx, 2, idx, A, INSERT_VALUES));
245
246 /*
247 Restore vector
248 */
249 PetscCall(VecRestoreArrayRead(x, &xx));
250
251 /*
252 Assemble matrix
253 */
254 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
255 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
256 PetscFunctionReturn(PETSC_SUCCESS);
257 }
258
FormFunction2(SNES snes,Vec x,Vec f,void * dummy)259 static PetscErrorCode FormFunction2(SNES snes, Vec x, Vec f, void *dummy)
260 {
261 const PetscScalar *xx;
262 PetscScalar *ff;
263
264 PetscFunctionBeginUser;
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 PetscCall(VecGetArrayRead(x, &xx));
273 PetscCall(VecGetArray(f, &ff));
274
275 /*
276 Compute function
277 */
278 ff[0] = PetscSinScalar(3.0 * xx[0]) + xx[0];
279 ff[1] = xx[1];
280
281 /*
282 Restore vectors
283 */
284 PetscCall(VecRestoreArrayRead(x, &xx));
285 PetscCall(VecRestoreArray(f, &ff));
286 PetscFunctionReturn(PETSC_SUCCESS);
287 }
288
FormJacobian2(SNES snes,Vec x,Mat jac,Mat B,void * dummy)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};
294
295 PetscFunctionBeginUser;
296 /*
297 Get pointer to vector data
298 */
299 PetscCall(VecGetArrayRead(x, &xx));
300
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 PetscCall(MatSetValues(jac, 2, idx, 2, idx, A, INSERT_VALUES));
311
312 /*
313 Restore vector
314 */
315 PetscCall(VecRestoreArrayRead(x, &xx));
316
317 /*
318 Assemble matrix
319 */
320 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
321 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
322 PetscFunctionReturn(PETSC_SUCCESS);
323 }
324
block_system(void)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;
335
336 Mat j11, j12, j21, j22;
337 Vec x1, x2, r1, r2;
338 Vec bv;
339 Vec bx[2];
340 Mat bA[2][2];
341
342 PetscFunctionBeginUser;
343 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n\n========================= Block system =========================\n\n"));
344
345 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
346
347 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
348 Create matrix and vector data structures; set corresponding routines
349 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
350
351 /*
352 Create sub vectors for solution and nonlinear function
353 */
354 PetscCall(VecCreateSeq(PETSC_COMM_SELF, 1, &x1));
355 PetscCall(VecDuplicate(x1, &r1));
356
357 PetscCall(VecCreateSeq(PETSC_COMM_SELF, 1, &x2));
358 PetscCall(VecDuplicate(x2, &r2));
359
360 /*
361 Create the block vectors
362 */
363 bx[0] = x1;
364 bx[1] = x2;
365 PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, NULL, bx, &x));
366 PetscCall(VecAssemblyBegin(x));
367 PetscCall(VecAssemblyEnd(x));
368 PetscCall(VecDestroy(&x1));
369 PetscCall(VecDestroy(&x2));
370
371 bx[0] = r1;
372 bx[1] = r2;
373 PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, NULL, bx, &r));
374 PetscCall(VecDestroy(&r1));
375 PetscCall(VecDestroy(&r2));
376 PetscCall(VecAssemblyBegin(r));
377 PetscCall(VecAssemblyEnd(r));
378
379 /*
380 Create sub Jacobian matrix data structure
381 */
382 PetscCall(MatCreate(PETSC_COMM_WORLD, &j11));
383 PetscCall(MatSetSizes(j11, 1, 1, 1, 1));
384 PetscCall(MatSetType(j11, MATSEQAIJ));
385 PetscCall(MatSetUp(j11));
386
387 PetscCall(MatCreate(PETSC_COMM_WORLD, &j12));
388 PetscCall(MatSetSizes(j12, 1, 1, 1, 1));
389 PetscCall(MatSetType(j12, MATSEQAIJ));
390 PetscCall(MatSetUp(j12));
391
392 PetscCall(MatCreate(PETSC_COMM_WORLD, &j21));
393 PetscCall(MatSetSizes(j21, 1, 1, 1, 1));
394 PetscCall(MatSetType(j21, MATSEQAIJ));
395 PetscCall(MatSetUp(j21));
396
397 PetscCall(MatCreate(PETSC_COMM_WORLD, &j22));
398 PetscCall(MatSetSizes(j22, PETSC_DECIDE, PETSC_DECIDE, 1, 1));
399 PetscCall(MatSetType(j22, MATSEQAIJ));
400 PetscCall(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;
408
409 PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, &bA[0][0], &J));
410 PetscCall(MatSetUp(J));
411 PetscCall(MatNestSetVecType(J, VECNEST));
412 PetscCall(MatDestroy(&j11));
413 PetscCall(MatDestroy(&j12));
414 PetscCall(MatDestroy(&j21));
415 PetscCall(MatDestroy(&j22));
416
417 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
418 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
419
420 PetscCall(PetscOptionsHasName(NULL, NULL, "-hard", &flg));
421 if (!flg) {
422 /*
423 Set function evaluation routine and vector.
424 */
425 PetscCall(SNESSetFunction(snes, r, FormFunction1_block, NULL));
426
427 /*
428 Set Jacobian matrix data structure and Jacobian evaluation routine
429 */
430 PetscCall(SNESSetJacobian(snes, J, J, FormJacobian1_block, NULL));
431 } else {
432 PetscCall(SNESSetFunction(snes, r, FormFunction2_block, NULL));
433 PetscCall(SNESSetJacobian(snes, J, J, FormJacobian2_block, NULL));
434 }
435
436 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
437 Customize nonlinear solver; set runtime options
438 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
439
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 PetscCall(SNESGetKSP(snes, &ksp));
446 PetscCall(KSPGetPC(ksp, &pc));
447 PetscCall(PCSetType(pc, PCNONE));
448 PetscCall(KSPSetTolerances(ksp, 1.e-4, PETSC_CURRENT, PETSC_CURRENT, 20));
449
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 PetscCall(SNESSetFromOptions(snes));
458
459 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
460 Evaluate initial guess; then solve nonlinear system
461 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
462 if (!flg) {
463 PetscCall(VecSet(x, pfive));
464 } else {
465 Vec *vecs;
466 PetscCall(VecNestGetSubVecs(x, NULL, &vecs));
467 bv = vecs[0];
468 /* PetscCall(VecBlockGetSubVec(x, 0, &bv)); */
469 PetscCall(VecSetValue(bv, 0, 2.0, INSERT_VALUES)); /* xx[0] = 2.0; */
470 PetscCall(VecAssemblyBegin(bv));
471 PetscCall(VecAssemblyEnd(bv));
472
473 /* PetscCall(VecBlockGetSubVec(x, 1, &bv)); */
474 bv = vecs[1];
475 PetscCall(VecSetValue(bv, 0, 3.0, INSERT_VALUES)); /* xx[1] = 3.0; */
476 PetscCall(VecAssemblyBegin(bv));
477 PetscCall(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 PetscCall(SNESSolve(snes, NULL, x));
486 PetscCall(SNESGetIterationNumber(snes, &its));
487 if (flg) {
488 Vec f;
489 PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
490 PetscCall(SNESGetFunction(snes, &f, 0, 0));
491 PetscCall(VecView(r, PETSC_VIEWER_STDOUT_WORLD));
492 }
493
494 PetscCall(PetscPrintf(PETSC_COMM_SELF, "number of SNES iterations = %" PetscInt_FMT "\n\n", its));
495
496 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
497 Free work space. All PETSc objects should be destroyed when they
498 are no longer needed.
499 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
500 PetscCall(VecDestroy(&x));
501 PetscCall(VecDestroy(&r));
502 PetscCall(MatDestroy(&J));
503 PetscCall(SNESDestroy(&snes));
504 PetscFunctionReturn(PETSC_SUCCESS);
505 }
506
FormFunction1_block(SNES snes,Vec x,Vec f,void * dummy)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;
513
514 PetscFunctionBeginUser;
515 /* get blocks for function */
516 PetscCall(VecNestGetSubVecs(f, &nb, &ff));
517 f1 = ff[0];
518 f2 = ff[1];
519
520 /* get blocks for solution */
521 PetscCall(VecNestGetSubVecs(x, &nb, &xx));
522 x1 = xx[0];
523 x2 = xx[1];
524
525 /* get solution values */
526 index = 0;
527 PetscCall(VecGetValues(x1, 1, &index, &xx_0));
528 PetscCall(VecGetValues(x2, 1, &index, &xx_1));
529
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;
533
534 /* set function values */
535 PetscCall(VecSetValue(f1, index, ff_0, INSERT_VALUES));
536
537 PetscCall(VecSetValue(f2, index, ff_1, INSERT_VALUES));
538
539 PetscCall(VecAssemblyBegin(f));
540 PetscCall(VecAssemblyEnd(f));
541 PetscFunctionReturn(PETSC_SUCCESS);
542 }
543
FormJacobian1_block(SNES snes,Vec x,Mat jac,Mat B,void * dummy)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;
552
553 PetscFunctionBeginUser;
554 /* get blocks for solution */
555 PetscCall(VecNestGetSubVecs(x, &nb, &xx));
556 x1 = xx[0];
557 x2 = xx[1];
558
559 /* get solution values */
560 index = 0;
561 PetscCall(VecGetValues(x1, 1, &index, &xx_0));
562 PetscCall(VecGetValues(x2, 1, &index, &xx_1));
563
564 /* get block matrices */
565 PetscCall(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];
570
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;
576
577 /* set jacobian values */
578 PetscCall(MatSetValue(j11, 0, 0, A_00, INSERT_VALUES));
579 PetscCall(MatSetValue(j12, 0, 0, A_01, INSERT_VALUES));
580 PetscCall(MatSetValue(j21, 0, 0, A_10, INSERT_VALUES));
581 PetscCall(MatSetValue(j22, 0, 0, A_11, INSERT_VALUES));
582
583 /* Assemble sub matrix */
584 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
585 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
586 PetscFunctionReturn(PETSC_SUCCESS);
587 }
588
FormFunction2_block(SNES snes,Vec x,Vec f,void * dummy)589 static PetscErrorCode FormFunction2_block(SNES snes, Vec x, Vec f, void *dummy)
590 {
591 PetscScalar *ff;
592 const PetscScalar *xx;
593
594 PetscFunctionBeginUser;
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 PetscCall(VecGetArrayRead(x, &xx));
603 PetscCall(VecGetArray(f, &ff));
604
605 /*
606 Compute function
607 */
608 ff[0] = PetscSinScalar(3.0 * xx[0]) + xx[0];
609 ff[1] = xx[1];
610
611 /*
612 Restore vectors
613 */
614 PetscCall(VecRestoreArrayRead(x, &xx));
615 PetscCall(VecRestoreArray(f, &ff));
616 PetscFunctionReturn(PETSC_SUCCESS);
617 }
618
FormJacobian2_block(SNES snes,Vec x,Mat jac,Mat B,void * dummy)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};
624
625 PetscFunctionBeginUser;
626 /*
627 Get pointer to vector data
628 */
629 PetscCall(VecGetArrayRead(x, &xx));
630
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 PetscCall(MatSetValues(jac, 2, idx, 2, idx, A, INSERT_VALUES));
641
642 /*
643 Restore vector
644 */
645 PetscCall(VecRestoreArrayRead(x, &xx));
646
647 /*
648 Assemble matrix
649 */
650 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
651 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
652 PetscFunctionReturn(PETSC_SUCCESS);
653 }
654
main(int argc,char ** argv)655 int main(int argc, char **argv)
656 {
657 PetscMPIInt size;
658
659 PetscFunctionBeginUser;
660 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
661 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
662 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
663 PetscCall(assembled_system());
664 PetscCall(block_system());
665 PetscCall(PetscFinalize());
666 return 0;
667 }
668
669 /*TEST
670
671 test:
672 args: -snes_monitor_short
673 requires: !single
674
675 TEST*/
676