Actual source code: ex2.c


  2: static char help[] = "Newton method to solve u'' + u^{2} = f, sequentially.\n\
  3: This example employs a user-defined monitoring routine.\n\n";

  5: /*
  6:    Include "petscdraw.h" so that we can use PETSc drawing routines.
  7:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
  8:    file automatically includes:
  9:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 10:      petscmat.h - matrices
 11:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 12:      petscviewer.h - viewers               petscpc.h  - preconditioners
 13:      petscksp.h   - linear solvers
 14: */

 16: #include <petscsnes.h>

 18: /*
 19:    User-defined routines
 20: */
 21: extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
 22: extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
 23: extern PetscErrorCode FormInitialGuess(Vec);
 24: extern PetscErrorCode Monitor(SNES, PetscInt, PetscReal, void *);

 26: /*
 27:    User-defined context for monitoring
 28: */
 29: typedef struct {
 30:   PetscViewer viewer;
 31: } MonitorCtx;

 33: int main(int argc, char **argv)
 34: {
 35:   SNES        snes;       /* SNES context */
 36:   Vec         x, r, F, U; /* vectors */
 37:   Mat         J;          /* Jacobian matrix */
 38:   MonitorCtx  monP;       /* monitoring context */
 39:   PetscInt    its, n = 5, i, maxit, maxf;
 40:   PetscMPIInt size;
 41:   PetscScalar h, xp, v, none = -1.0;
 42:   PetscReal   abstol, rtol, stol, norm;

 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_SELF, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
 48:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
 49:   h = 1.0 / (n - 1);

 51:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 52:      Create nonlinear solver context
 53:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 55:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));

 57:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 58:      Create vector data structures; set function evaluation routine
 59:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 60:   /*
 61:      Note that we form 1 vector from scratch and then duplicate as needed.
 62:   */
 63:   PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
 64:   PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
 65:   PetscCall(VecSetFromOptions(x));
 66:   PetscCall(VecDuplicate(x, &r));
 67:   PetscCall(VecDuplicate(x, &F));
 68:   PetscCall(VecDuplicate(x, &U));

 70:   /*
 71:      Set function evaluation routine and vector
 72:   */
 73:   PetscCall(SNESSetFunction(snes, r, FormFunction, (void *)F));

 75:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 76:      Create matrix data structure; set Jacobian evaluation routine
 77:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 79:   PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
 80:   PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, n, n));
 81:   PetscCall(MatSetFromOptions(J));
 82:   PetscCall(MatSeqAIJSetPreallocation(J, 3, NULL));

 84:   /*
 85:      Set Jacobian matrix data structure and default Jacobian evaluation
 86:      routine. User can override with:
 87:      -snes_fd : default finite differencing approximation of Jacobian
 88:      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
 89:                 (unless user explicitly sets preconditioner)
 90:      -snes_mf_operator : form preconditioning matrix as set by the user,
 91:                          but use matrix-free approx for Jacobian-vector
 92:                          products within Newton-Krylov method
 93:   */

 95:   PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, NULL));

 97:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 98:      Customize nonlinear solver; set runtime options
 99:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

101:   /*
102:      Set an optional user-defined monitoring routine
103:   */
104:   PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, 0, 0, 0, 400, 400, &monP.viewer));
105:   PetscCall(SNESMonitorSet(snes, Monitor, &monP, 0));

107:   /*
108:      Set names for some vectors to facilitate monitoring (optional)
109:   */
110:   PetscCall(PetscObjectSetName((PetscObject)x, "Approximate Solution"));
111:   PetscCall(PetscObjectSetName((PetscObject)U, "Exact Solution"));

113:   /*
114:      Set SNES/KSP/KSP/PC runtime options, e.g.,
115:          -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
116:   */
117:   PetscCall(SNESSetFromOptions(snes));

119:   /*
120:      Print parameters used for convergence testing (optional) ... just
121:      to demonstrate this routine; this information is also printed with
122:      the option -snes_view
123:   */
124:   PetscCall(SNESGetTolerances(snes, &abstol, &rtol, &stol, &maxit, &maxf));
125:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "atol=%g, rtol=%g, stol=%g, maxit=%" PetscInt_FMT ", maxf=%" PetscInt_FMT "\n", (double)abstol, (double)rtol, (double)stol, maxit, maxf));

127:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128:      Initialize application:
129:      Store right-hand-side of PDE and exact solution
130:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

132:   xp = 0.0;
133:   for (i = 0; i < n; i++) {
134:     v = 6.0 * xp + PetscPowScalar(xp + 1.e-12, 6.0); /* +1.e-12 is to prevent 0^6 */
135:     PetscCall(VecSetValues(F, 1, &i, &v, INSERT_VALUES));
136:     v = xp * xp * xp;
137:     PetscCall(VecSetValues(U, 1, &i, &v, INSERT_VALUES));
138:     xp += h;
139:   }

141:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
142:      Evaluate initial guess; then solve nonlinear system
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:   */
150:   PetscCall(FormInitialGuess(x));
151:   PetscCall(SNESSolve(snes, NULL, x));
152:   PetscCall(SNESGetIterationNumber(snes, &its));
153:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "number of SNES iterations = %" PetscInt_FMT "\n\n", its));

155:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156:      Check solution and clean up
157:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

159:   /*
160:      Check the error
161:   */
162:   PetscCall(VecAXPY(x, none, U));
163:   PetscCall(VecNorm(x, NORM_2, &norm));
164:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its));

166:   /*
167:      Free work space.  All PETSc objects should be destroyed when they
168:      are no longer needed.
169:   */
170:   PetscCall(VecDestroy(&x));
171:   PetscCall(VecDestroy(&r));
172:   PetscCall(VecDestroy(&U));
173:   PetscCall(VecDestroy(&F));
174:   PetscCall(MatDestroy(&J));
175:   PetscCall(SNESDestroy(&snes));
176:   PetscCall(PetscViewerDestroy(&monP.viewer));
177:   PetscCall(PetscFinalize());
178:   return 0;
179: }
180: /* ------------------------------------------------------------------- */
181: /*
182:    FormInitialGuess - Computes initial guess.

184:    Input/Output Parameter:
185: .  x - the solution vector
186: */
187: PetscErrorCode FormInitialGuess(Vec x)
188: {
189:   PetscScalar pfive = .50;
190:   PetscFunctionBeginUser;
191:   PetscCall(VecSet(x, pfive));
192:   PetscFunctionReturn(PETSC_SUCCESS);
193: }
194: /* ------------------------------------------------------------------- */
195: /*
196:    FormFunction - Evaluates nonlinear function, F(x).

198:    Input Parameters:
199: .  snes - the SNES context
200: .  x - input vector
201: .  ctx - optional user-defined context, as set by SNESSetFunction()

203:    Output Parameter:
204: .  f - function vector

206:    Note:
207:    The user-defined context can contain any application-specific data
208:    needed for the function evaluation (such as various parameters, work
209:    vectors, and grid information).  In this program the context is just
210:    a vector containing the right-hand-side of the discretized PDE.
211:  */

213: PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void *ctx)
214: {
215:   Vec                g = (Vec)ctx;
216:   const PetscScalar *xx, *gg;
217:   PetscScalar       *ff, d;
218:   PetscInt           i, n;

220:   PetscFunctionBeginUser;
221:   /*
222:      Get pointers to vector data.
223:        - For default PETSc vectors, VecGetArray() returns a pointer to
224:          the data array.  Otherwise, the routine is implementation dependent.
225:        - You MUST call VecRestoreArray() when you no longer need access to
226:          the array.
227:   */
228:   PetscCall(VecGetArrayRead(x, &xx));
229:   PetscCall(VecGetArray(f, &ff));
230:   PetscCall(VecGetArrayRead(g, &gg));

232:   /*
233:      Compute function
234:   */
235:   PetscCall(VecGetSize(x, &n));
236:   d     = (PetscReal)(n - 1);
237:   d     = d * d;
238:   ff[0] = xx[0];
239:   for (i = 1; i < n - 1; i++) ff[i] = d * (xx[i - 1] - 2.0 * xx[i] + xx[i + 1]) + xx[i] * xx[i] - gg[i];
240:   ff[n - 1] = xx[n - 1] - 1.0;

242:   /*
243:      Restore vectors
244:   */
245:   PetscCall(VecRestoreArrayRead(x, &xx));
246:   PetscCall(VecRestoreArray(f, &ff));
247:   PetscCall(VecRestoreArrayRead(g, &gg));
248:   PetscFunctionReturn(PETSC_SUCCESS);
249: }
250: /* ------------------------------------------------------------------- */
251: /*
252:    FormJacobian - Evaluates Jacobian matrix.

254:    Input Parameters:
255: .  snes - the SNES context
256: .  x - input vector
257: .  dummy - optional user-defined context (not used here)

259:    Output Parameters:
260: .  jac - Jacobian matrix
261: .  B - optionally different preconditioning matrix

263: */

265: PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
266: {
267:   const PetscScalar *xx;
268:   PetscScalar        A[3], d;
269:   PetscInt           i, n, j[3];

271:   PetscFunctionBeginUser;
272:   /*
273:      Get pointer to vector data
274:   */
275:   PetscCall(VecGetArrayRead(x, &xx));

277:   /*
278:      Compute Jacobian entries and insert into matrix.
279:       - Note that in this case we set all elements for a particular
280:         row at once.
281:   */
282:   PetscCall(VecGetSize(x, &n));
283:   d = (PetscReal)(n - 1);
284:   d = d * d;

286:   /*
287:      Interior grid points
288:   */
289:   for (i = 1; i < n - 1; i++) {
290:     j[0] = i - 1;
291:     j[1] = i;
292:     j[2] = i + 1;
293:     A[0] = A[2] = d;
294:     A[1]        = -2.0 * d + 2.0 * xx[i];
295:     PetscCall(MatSetValues(B, 1, &i, 3, j, A, INSERT_VALUES));
296:   }

298:   /*
299:      Boundary points
300:   */
301:   i    = 0;
302:   A[0] = 1.0;

304:   PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES));

306:   i    = n - 1;
307:   A[0] = 1.0;

309:   PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES));

311:   /*
312:      Restore vector
313:   */
314:   PetscCall(VecRestoreArrayRead(x, &xx));

316:   /*
317:      Assemble matrix
318:   */
319:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
320:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
321:   if (jac != B) {
322:     PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
323:     PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
324:   }
325:   PetscFunctionReturn(PETSC_SUCCESS);
326: }
327: /* ------------------------------------------------------------------- */
328: /*
329:    Monitor - User-defined monitoring routine that views the
330:    current iterate with an x-window plot.

332:    Input Parameters:
333:    snes - the SNES context
334:    its - iteration number
335:    norm - 2-norm function value (may be estimated)
336:    ctx - optional user-defined context for private data for the
337:          monitor routine, as set by SNESMonitorSet()

339:    Note:
340:    See the manpage for PetscViewerDrawOpen() for useful runtime options,
341:    such as -nox to deactivate all x-window output.
342:  */
343: PetscErrorCode Monitor(SNES snes, PetscInt its, PetscReal fnorm, void *ctx)
344: {
345:   MonitorCtx *monP = (MonitorCtx *)ctx;
346:   Vec         x;

348:   PetscFunctionBeginUser;
349:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "iter = %" PetscInt_FMT ", SNES Function norm %g\n", its, (double)fnorm));
350:   PetscCall(SNESGetSolution(snes, &x));
351:   PetscCall(VecView(x, monP->viewer));
352:   PetscFunctionReturn(PETSC_SUCCESS);
353: }

355: /*TEST

357:    test:
358:       args: -nox -snes_monitor_cancel -snes_monitor_short -snes_view -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always

360:    test:
361:       suffix: 2
362:       args: -nox -snes_monitor_cancel -snes_monitor_short -snes_type newtontr -snes_view
363:       requires: !single

365:    test:
366:       suffix: 3
367:       args: -nox -malloc no -options_left no -snes_monitor_cancel -snes_monitor_short -snes_view -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always

369:    test:
370:       suffix: 4
371:       args: -nox -snes_monitor_cancel -snes_monitor_short -snes_type newtontrdc -snes_view
372:       requires: !single

374: TEST*/