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:   PetscInitialize(&argc,&argv,(char*)0,help);
 45:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 47:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 48:   h    = 1.0/(n-1);

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

 54:   SNESCreate(PETSC_COMM_WORLD,&snes);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

180:    Input/Output Parameter:
181: .  x - the solution vector
182: */
183: PetscErrorCode FormInitialGuess(Vec x)
184: {
185:   PetscScalar    pfive = .50;
186:   VecSet(x,pfive);
187:   return 0;
188: }
189: /* ------------------------------------------------------------------- */
190: /*
191:    FormFunction - Evaluates nonlinear function, F(x).

193:    Input Parameters:
194: .  snes - the SNES context
195: .  x - input vector
196: .  ctx - optional user-defined context, as set by SNESSetFunction()

198:    Output Parameter:
199: .  f - function vector

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

208: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
209: {
210:   Vec               g = (Vec)ctx;
211:   const PetscScalar *xx,*gg;
212:   PetscScalar       *ff,d;
213:   PetscInt          i,n;

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

226:   /*
227:      Compute function
228:   */
229:   VecGetSize(x,&n);
230:   d     = (PetscReal)(n - 1); d = d*d;
231:   ff[0] = xx[0];
232:   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];
233:   ff[n-1] = xx[n-1] - 1.0;

235:   /*
236:      Restore vectors
237:   */
238:   VecRestoreArrayRead(x,&xx);
239:   VecRestoreArray(f,&ff);
240:   VecRestoreArrayRead(g,&gg);
241:   return 0;
242: }
243: /* ------------------------------------------------------------------- */
244: /*
245:    FormJacobian - Evaluates Jacobian matrix.

247:    Input Parameters:
248: .  snes - the SNES context
249: .  x - input vector
250: .  dummy - optional user-defined context (not used here)

252:    Output Parameters:
253: .  jac - Jacobian matrix
254: .  B - optionally different preconditioning matrix

256: */

258: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
259: {
260:   const PetscScalar *xx;
261:   PetscScalar       A[3],d;
262:   PetscInt          i,n,j[3];

264:   /*
265:      Get pointer to vector data
266:   */
267:   VecGetArrayRead(x,&xx);

269:   /*
270:      Compute Jacobian entries and insert into matrix.
271:       - Note that in this case we set all elements for a particular
272:         row at once.
273:   */
274:   VecGetSize(x,&n);
275:   d    = (PetscReal)(n - 1); d = d*d;

277:   /*
278:      Interior grid points
279:   */
280:   for (i=1; i<n-1; i++) {
281:     j[0] = i - 1; j[1] = i; j[2] = i + 1;
282:     A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
283:     MatSetValues(B,1,&i,3,j,A,INSERT_VALUES);
284:   }

286:   /*
287:      Boundary points
288:   */
289:   i = 0;   A[0] = 1.0;

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

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

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

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

302:   /*
303:      Assemble matrix
304:   */
305:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
306:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
307:   if (jac != B) {
308:     MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
309:     MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
310:   }
311:   return 0;
312: }
313: /* ------------------------------------------------------------------- */
314: /*
315:    Monitor - User-defined monitoring routine that views the
316:    current iterate with an x-window plot.

318:    Input Parameters:
319:    snes - the SNES context
320:    its - iteration number
321:    norm - 2-norm function value (may be estimated)
322:    ctx - optional user-defined context for private data for the
323:          monitor routine, as set by SNESMonitorSet()

325:    Note:
326:    See the manpage for PetscViewerDrawOpen() for useful runtime options,
327:    such as -nox to deactivate all x-window output.
328:  */
329: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
330: {
331:   MonitorCtx     *monP = (MonitorCtx*) ctx;
332:   Vec            x;

334:   PetscPrintf(PETSC_COMM_WORLD,"iter = %D, SNES Function norm %g\n",its,(double)fnorm);
335:   SNESGetSolution(snes,&x);
336:   VecView(x,monP->viewer);
337:   return 0;
338: }

340: /*TEST

342:    test:
343:       args: -nox -snes_monitor_cancel -snes_monitor_short -snes_view -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always

345:    test:
346:       suffix: 2
347:       args: -nox -snes_monitor_cancel -snes_monitor_short -snes_type newtontr -snes_view
348:       requires: !single

350:    test:
351:       suffix: 3
352:       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

354:    test:
355:       suffix: 4
356:       args: -nox -snes_monitor_cancel -snes_monitor_short -snes_type newtontrdc -snes_view
357:       requires: !single

359: TEST*/