Actual source code: ex1.c


  2: static char help[] = "Newton's method for a two-variable system, sequential.\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:      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: /*F
 14: This examples solves either
 15: \begin{equation}
 16:   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}
 17: \end{equation}
 18: or if the {\tt -hard} options is given
 19: \begin{equation}
 20:   F\genfrac{(}{)}{0pt}{}{x_0}{x_1} = \genfrac{(}{)}{0pt}{}{\sin(3 x_0) + x_0}{x_1}
 21: \end{equation}
 22: F*/
 23: #include <petscsnes.h>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

208:   PetscFunctionBeginUser;
209:   /*
210:      Get pointer to vector data
211:   */
212:   PetscCall(VecGetArrayRead(x, &xx));

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

225:   /*
226:      Restore vector
227:   */
228:   PetscCall(VecRestoreArrayRead(x, &xx));

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

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

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

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

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

279:   PetscFunctionBeginUser;
280:   /*
281:      Get pointer to vector data
282:   */
283:   PetscCall(VecGetArrayRead(x, &xx));

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

296:   /*
297:      Restore vector
298:   */
299:   PetscCall(VecRestoreArrayRead(x, &xx));

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

313: /*TEST

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

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

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

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

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

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

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

355: TEST*/