Actual source code: ex1.c

  1: static char help[] = "Newton's method for a two-variable system, sequential.\n\n";

3: /*
4:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
5:    file automatically includes:
6:      petscsys.h       - base PETSc routines   petscvec.h - vectors
7:      petscmat.h - matrices
8:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
9:      petscviewer.h - viewers               petscpc.h  - preconditioners
10:      petscksp.h   - linear solvers
11: */
12: /*F
13: This examples solves either
14:
15:   F\genfrac{(}{)}{0pt}{}{x_0}{x_1} = \genfrac{(}{)}{0pt}{}{x^2_0 + x_0 x_1 - 3}{x_0 x_1 + x^2_1 - 6}
16:
17: or if the {\tt -hard} options is given
18:
19:   F\genfrac{(}{)}{0pt}{}{x_0}{x_1} = \genfrac{(}{)}{0pt}{}{\sin(3 x_0) + x_0}{x_1}
20:
21: F*/
22: #include <petscsnes.h>

24: /*
25:    User-defined routines
26: */
27: extern PetscErrorCode FormJacobian1(SNES, Vec, Mat, Mat, void *);
28: extern PetscErrorCode FormFunction1(SNES, Vec, Vec, void *);
29: extern PetscErrorCode FormJacobian2(SNES, Vec, Mat, Mat, void *);
30: extern PetscErrorCode FormFunction2(SNES, Vec, Vec, void *);

32: int main(int argc, char **argv)
33: {
34:   SNES        snes; /* nonlinear solver context */
35:   KSP         ksp;  /* linear solver context */
36:   PC          pc;   /* preconditioner context */
37:   Vec         x, r; /* solution, residual vectors */
38:   Mat         J;    /* Jacobian matrix */
39:   PetscMPIInt size;
40:   PetscScalar pfive = .5, *xx;
41:   PetscBool   flg;

43:   PetscFunctionBeginUser;
44:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
45:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
46:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Example is only for sequential runs");

48:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
49:      Create nonlinear solver context
50:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
51:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
52:   PetscCall(SNESSetType(snes, SNESNEWTONLS));
53:   PetscCall(SNESSetOptionsPrefix(snes, "mysolver_"));

55:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
56:      Create matrix and vector data structures; set corresponding routines
57:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
58:   /*
59:      Create vectors for solution and nonlinear function
60:   */
61:   PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
62:   PetscCall(VecSetSizes(x, PETSC_DECIDE, 2));
63:   PetscCall(VecSetFromOptions(x));
64:   PetscCall(VecDuplicate(x, &r));

66:   /*
67:      Create Jacobian matrix data structure
68:   */
69:   PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
70:   PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, 2, 2));
71:   PetscCall(MatSetFromOptions(J));
72:   PetscCall(MatSetUp(J));

74:   PetscCall(PetscOptionsHasName(NULL, NULL, "-hard", &flg));
75:   if (!flg) {
76:     /*
77:      Set function evaluation routine and vector.
78:     */
79:     PetscCall(SNESSetFunction(snes, r, FormFunction1, NULL));

81:     /*
82:      Set Jacobian matrix data structure and Jacobian evaluation routine
83:     */
84:     PetscCall(SNESSetJacobian(snes, J, J, FormJacobian1, NULL));
85:   } else {
86:     PetscCall(SNESSetFunction(snes, r, FormFunction2, NULL));
87:     PetscCall(SNESSetJacobian(snes, J, J, FormJacobian2, NULL));
88:   }

90:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91:      Customize nonlinear solver; set runtime options
92:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
93:   /*
94:      Set linear solver defaults for this problem. By extracting the
95:      KSP and PC contexts from the SNES context, we can then
96:      directly call any KSP and PC routines to set various options.
97:   */
98:   PetscCall(SNESGetKSP(snes, &ksp));
99:   PetscCall(KSPGetPC(ksp, &pc));
100:   PetscCall(PCSetType(pc, PCNONE));
101:   PetscCall(KSPSetTolerances(ksp, 1.e-4, PETSC_DEFAULT, PETSC_DEFAULT, 20));

103:   /*
104:      Set SNES/KSP/KSP/PC runtime options, e.g.,
105:          -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
106:      These options will override those specified above as long as
107:      SNESSetFromOptions() is called _after_ any other customization
108:      routines.
109:   */
110:   PetscCall(SNESSetFromOptions(snes));

112:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113:      Evaluate initial guess; then solve nonlinear system
114:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
115:   if (!flg) {
116:     PetscCall(VecSet(x, pfive));
117:   } else {
118:     PetscCall(VecGetArray(x, &xx));
119:     xx[0] = 2.0;
120:     xx[1] = 3.0;
121:     PetscCall(VecRestoreArray(x, &xx));
122:   }
123:   /*
124:      Note: The user should initialize the vector, x, with the initial guess
125:      for the nonlinear solver prior to calling SNESSolve().  In particular,
126:      to employ an initial guess of zero, the user should explicitly set
127:      this vector to zero by calling VecSet().
128:   */

130:   PetscCall(SNESSolve(snes, NULL, x));
131:   if (flg) {
132:     Vec f;
133:     PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
134:     PetscCall(SNESGetFunction(snes, &f, 0, 0));
135:     PetscCall(VecView(r, PETSC_VIEWER_STDOUT_WORLD));
136:   }

138:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
139:      Free work space.  All PETSc objects should be destroyed when they
140:      are no longer needed.
141:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

143:   PetscCall(VecDestroy(&x));
144:   PetscCall(VecDestroy(&r));
145:   PetscCall(MatDestroy(&J));
146:   PetscCall(SNESDestroy(&snes));
147:   PetscCall(PetscFinalize());
148:   return 0;
149: }
150: /* ------------------------------------------------------------------- */
151: /*
152:    FormFunction1 - Evaluates nonlinear function, F(x).

154:    Input Parameters:
155: .  snes - the SNES context
156: .  x    - input vector
157: .  ctx  - optional user-defined context

159:    Output Parameter:
160: .  f - function vector
161:  */
162: PetscErrorCode FormFunction1(SNES snes, Vec x, Vec f, void *ctx)
163: {
164:   const PetscScalar *xx;
165:   PetscScalar       *ff;

167:   PetscFunctionBeginUser;
168:   /*
169:    Get pointers to vector data.
170:       - For default PETSc vectors, VecGetArray() returns a pointer to
171:         the data array.  Otherwise, the routine is implementation dependent.
172:       - You MUST call VecRestoreArray() when you no longer need access to
173:         the array.
174:    */
176:   PetscCall(VecGetArray(f, &ff));

178:   /* Compute function */
179:   ff[0] = xx[0] * xx[0] + xx[0] * xx[1] - 3.0;
180:   ff[1] = xx[0] * xx[1] + xx[1] * xx[1] - 6.0;

182:   /* Restore vectors */
184:   PetscCall(VecRestoreArray(f, &ff));
185:   PetscFunctionReturn(PETSC_SUCCESS);
186: }
187: /* ------------------------------------------------------------------- */
188: /*
189:    FormJacobian1 - Evaluates Jacobian matrix.

191:    Input Parameters:
192: .  snes - the SNES context
193: .  x - input vector
194: .  dummy - optional user-defined context (not used here)

196:    Output Parameters:
197: .  jac - Jacobian matrix
198: .  B - optionally different preconditioning matrix
199: .  flag - flag indicating matrix structure
200: */
201: PetscErrorCode FormJacobian1(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
202: {
203:   const PetscScalar *xx;
204:   PetscScalar        A[4];
205:   PetscInt           idx[2] = {0, 1};

207:   PetscFunctionBeginUser;
208:   /*
209:      Get pointer to vector data
210:   */

213:   /*
214:      Compute Jacobian entries and insert into matrix.
215:       - Since this is such a small problem, we set all entries for
216:         the matrix at once.
217:   */
218:   A[0] = 2.0 * xx[0] + xx[1];
219:   A[1] = xx[0];
220:   A[2] = xx[1];
221:   A[3] = xx[0] + 2.0 * xx[1];
222:   PetscCall(MatSetValues(B, 2, idx, 2, idx, A, INSERT_VALUES));

224:   /*
225:      Restore vector
226:   */

229:   /*
230:      Assemble matrix
231:   */
232:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
233:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
234:   if (jac != B) {
235:     PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
236:     PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
237:   }
238:   PetscFunctionReturn(PETSC_SUCCESS);
239: }

241: /* ------------------------------------------------------------------- */
242: PetscErrorCode FormFunction2(SNES snes, Vec x, Vec f, void *dummy)
243: {
244:   const PetscScalar *xx;
245:   PetscScalar       *ff;

247:   PetscFunctionBeginUser;
248:   /*
249:      Get pointers to vector data.
250:        - For default PETSc vectors, VecGetArray() returns a pointer to
251:          the data array.  Otherwise, the routine is implementation dependent.
252:        - You MUST call VecRestoreArray() when you no longer need access to
253:          the array.
254:   */
256:   PetscCall(VecGetArray(f, &ff));

258:   /*
259:      Compute function
260:   */
261:   ff[0] = PetscSinScalar(3.0 * xx[0]) + xx[0];
262:   ff[1] = xx[1];

264:   /*
265:      Restore vectors
266:   */
268:   PetscCall(VecRestoreArray(f, &ff));
269:   PetscFunctionReturn(PETSC_SUCCESS);
270: }
271: /* ------------------------------------------------------------------- */
272: PetscErrorCode FormJacobian2(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
273: {
274:   const PetscScalar *xx;
275:   PetscScalar        A[4];
276:   PetscInt           idx[2] = {0, 1};

278:   PetscFunctionBeginUser;
279:   /*
280:      Get pointer to vector data
281:   */

284:   /*
285:      Compute Jacobian entries and insert into matrix.
286:       - Since this is such a small problem, we set all entries for
287:         the matrix at once.
288:   */
289:   A[0] = 3.0 * PetscCosScalar(3.0 * xx[0]) + 1.0;
290:   A[1] = 0.0;
291:   A[2] = 0.0;
292:   A[3] = 1.0;
293:   PetscCall(MatSetValues(B, 2, idx, 2, idx, A, INSERT_VALUES));

295:   /*
296:      Restore vector
297:   */

300:   /*
301:      Assemble matrix
302:   */
303:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
304:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
305:   if (jac != B) {
306:     PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
307:     PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
308:   }
309:   PetscFunctionReturn(PETSC_SUCCESS);
310: }

312: /*TEST

314:    test:
315:       args: -prefix_push mysolver_ -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short -prefix_pop
316:       requires: !single

318:    # test harness puts {{ }} options always at the end, need to specify the prefix explicitly
319:    test:
320:       suffix: 2
321:       requires: !single
322:       args: -prefix_push mysolver_ -snes_monitor_short -prefix_pop -mysolver_snes_ksp_ew {{0 1}}
323:       output_file: output/ex1_1.out

325:    test:
326:       suffix: 3
327:       args: -prefix_push mysolver_ -ksp_view_solution ascii:ex1_2_sol.tmp:ascii_matlab -snes_monitor_short -prefix_pop
328:       requires: !single
329:       output_file: output/ex1_1.out

331:    test:
332:       suffix: 4
333:       args: -prefix_push mysolver_ -ksp_view_solution ascii:ex1_2_sol.tmp::append -snes_monitor_short -prefix_pop
334:       requires: !single
335:       output_file: output/ex1_1.out

337:    test:
338:       suffix: 5
339:       args: -prefix_push mysolver_ -ksp_view_solution ascii:ex1_2_sol.tmp:ascii_matlab:append -snes_monitor_short -prefix_pop
340:       requires: !single
341:       output_file: output/ex1_1.out

343:    test:
344:       suffix: 6
345:       args: -prefix_push mysolver_ -ksp_view_solution ascii:ex1_2_sol.tmp:default:append -snes_monitor_short -prefix_pop
346:       requires: !single
347:       output_file: output/ex1_1.out

349:    test:
350:       suffix: X
351:       args: -prefix_push mysolver_ -ksp_monitor_short -ksp_type gmres -ksp_gmres_krylov_monitor -snes_monitor_short -snes_rtol 1.e-4 -prefix_pop
352:       requires: !single x

354: TEST*/