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: /*T
  6:    Concepts: SNES^basic uniprocessor example
  7:    Concepts: SNES^setting a user-defined monitoring routine
  8:    Processors: 1
  9: T*/

 11: /*
 12:    Include "petscdraw.h" so that we can use PETSc drawing routines.
 13:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
 14:    file automatically includes:
 15:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 16:      petscmat.h - matrices
 17:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 18:      petscviewer.h - viewers               petscpc.h  - preconditioners
 19:      petscksp.h   - linear solvers
 20: */

 22: #include <petscsnes.h>

 24: /*
 25:    User-defined routines
 26: */
 27: extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
 28: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
 29: extern PetscErrorCode FormInitialGuess(Vec);
 30: extern PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*);

 32: /*
 33:    User-defined context for monitoring
 34: */
 35: typedef struct {
 36:   PetscViewer viewer;
 37: } MonitorCtx;

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

 50:   PetscInitialize(&argc,&argv,(char*)0,help);
 51:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 53:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 54:   h    = 1.0/(n-1);

 56:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 57:      Create nonlinear solver context
 58:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 60:   SNESCreate(PETSC_COMM_WORLD,&snes);

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

 75:   /*
 76:      Set function evaluation routine and vector
 77:   */
 78:   SNESSetFunction(snes,r,FormFunction,(void*)F);

 80:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 81:      Create matrix data structure; set Jacobian evaluation routine
 82:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 84:   MatCreate(PETSC_COMM_WORLD,&J);
 85:   MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);
 86:   MatSetFromOptions(J);
 87:   MatSeqAIJSetPreallocation(J,3,NULL);

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

100:   SNESSetJacobian(snes,J,J,FormJacobian,NULL);

102:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103:      Customize nonlinear solver; set runtime options
104:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

106:   /*
107:      Set an optional user-defined monitoring routine
108:   */
109:   PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,400,400,&monP.viewer);
110:   SNESMonitorSet(snes,Monitor,&monP,0);

112:   /*
113:      Set names for some vectors to facilitate monitoring (optional)
114:   */
115:   PetscObjectSetName((PetscObject)x,"Approximate Solution");
116:   PetscObjectSetName((PetscObject)U,"Exact Solution");

118:   /*
119:      Set SNES/KSP/KSP/PC runtime options, e.g.,
120:          -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
121:   */
122:   SNESSetFromOptions(snes);

124:   /*
125:      Print parameters used for convergence testing (optional) ... just
126:      to demonstrate this routine; this information is also printed with
127:      the option -snes_view
128:   */
129:   SNESGetTolerances(snes,&abstol,&rtol,&stol,&maxit,&maxf);
130:   PetscPrintf(PETSC_COMM_WORLD,"atol=%g, rtol=%g, stol=%g, maxit=%D, maxf=%D\n",(double)abstol,(double)rtol,(double)stol,maxit,maxf);

132:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
133:      Initialize application:
134:      Store right-hand-side of PDE and exact solution
135:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

137:   xp = 0.0;
138:   for (i=0; i<n; i++) {
139:     v    = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
140:     VecSetValues(F,1,&i,&v,INSERT_VALUES);
141:     v    = xp*xp*xp;
142:     VecSetValues(U,1,&i,&v,INSERT_VALUES);
143:     xp  += h;
144:   }

146:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147:      Evaluate initial guess; then solve nonlinear system
148:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
149:   /*
150:      Note: The user should initialize the vector, x, with the initial guess
151:      for the nonlinear solver prior to calling SNESSolve().  In particular,
152:      to employ an initial guess of zero, the user should explicitly set
153:      this vector to zero by calling VecSet().
154:   */
155:   FormInitialGuess(x);
156:   SNESSolve(snes,NULL,x);
157:   SNESGetIterationNumber(snes,&its);
158:   PetscPrintf(PETSC_COMM_WORLD,"number of SNES iterations = %D\n\n",its);

160:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
161:      Check solution and clean up
162:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

164:   /*
165:      Check the error
166:   */
167:   VecAXPY(x,none,U);
168:   VecNorm(x,NORM_2,&norm);
169:   PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);

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

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

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

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

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

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

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:   VecGetArrayRead(x,&xx);
229:   VecGetArray(f,&ff);
230:   VecGetArrayRead(g,&gg);

232:   /*
233:      Compute function
234:   */
235:   VecGetSize(x,&n);
236:   d     = (PetscReal)(n - 1); d = d*d;
237:   ff[0] = xx[0];
238:   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];
239:   ff[n-1] = xx[n-1] - 1.0;

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

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

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

262: */

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

270:   /*
271:      Get pointer to vector data
272:   */
273:   VecGetArrayRead(x,&xx);

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

283:   /*
284:      Interior grid points
285:   */
286:   for (i=1; i<n-1; i++) {
287:     j[0] = i - 1; j[1] = i; j[2] = i + 1;
288:     A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
289:     MatSetValues(B,1,&i,3,j,A,INSERT_VALUES);
290:   }

292:   /*
293:      Boundary points
294:   */
295:   i = 0;   A[0] = 1.0;

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

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

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

303:   /*
304:      Restore vector
305:   */
306:   VecRestoreArrayRead(x,&xx);

308:   /*
309:      Assemble matrix
310:   */
311:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
312:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
313:   if (jac != B) {
314:     MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
315:     MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
316:   }
317:   return 0;
318: }
319: /* ------------------------------------------------------------------- */
320: /*
321:    Monitor - User-defined monitoring routine that views the
322:    current iterate with an x-window plot.

324:    Input Parameters:
325:    snes - the SNES context
326:    its - iteration number
327:    norm - 2-norm function value (may be estimated)
328:    ctx - optional user-defined context for private data for the
329:          monitor routine, as set by SNESMonitorSet()

331:    Note:
332:    See the manpage for PetscViewerDrawOpen() for useful runtime options,
333:    such as -nox to deactivate all x-window output.
334:  */
335: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
336: {
337:   MonitorCtx     *monP = (MonitorCtx*) ctx;
338:   Vec            x;

340:   PetscPrintf(PETSC_COMM_WORLD,"iter = %D, SNES Function norm %g\n",its,(double)fnorm);
341:   SNESGetSolution(snes,&x);
342:   VecView(x,monP->viewer);
343:   return 0;
344: }

346: /*TEST

348:    test:
349:       args: -nox -snes_monitor_cancel -snes_monitor_short -snes_view -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always

351:    test:
352:       suffix: 2
353:       args: -nox -snes_monitor_cancel -snes_monitor_short -snes_type newtontr -snes_view
354:       requires: !single

356:    test:
357:       suffix: 3
358:       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

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

365: TEST*/