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:      petsc.h       - base PETSc routines   petscvec.h - vectors
 16:      petscsys.h    - system routines       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: */

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

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

 53:   PetscInitialize(&argc,&argv,(char *)0,help);
 54:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 55:   if (size != 1) SETERRQ(1,"This is a uniprocessor example only!");
 56:   PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
 57:   h = 1.0/(n-1);

 59:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 60:      Create nonlinear solver context
 61:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 63:   SNESCreate(PETSC_COMM_WORLD,&snes);

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

 78:   /* 
 79:      Set function evaluation routine and vector
 80:   */
 81:   SNESSetFunction(snes,r,FormFunction,(void*)F);


 84:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 85:      Create matrix data structure; set Jacobian evaluation routine
 86:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 88:   MatCreate(PETSC_COMM_WORLD,&J);
 89:   MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);
 90:   MatSetFromOptions(J);

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

103:   SNESSetJacobian(snes,J,J,FormJacobian,PETSC_NULL);

105:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106:      Customize nonlinear solver; set runtime options
107:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

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

135:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136:      Initialize application:
137:      Store right-hand-side of PDE and exact solution
138:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

163:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164:      Check solution and clean up
165:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

167:   /* 
168:      Check the error
169:   */
170:   VecAXPY(x,none,U);
171:   VecNorm(x,NORM_2,&norm);
172:   PetscPrintf(PETSC_COMM_WORLD,"Norm of error %A, Iterations %D\n",norm,its);


175:   /*
176:      Free work space.  All PETSc objects should be destroyed when they
177:      are no longer needed.
178:   */
179:   VecDestroy(x);  VecDestroy(r);
180:   VecDestroy(U);  VecDestroy(F);
181:   MatDestroy(J);  SNESDestroy(snes);
182:   PetscViewerDestroy(monP.viewer);
183:   PetscFinalize();

185:   return 0;
186: }
187: /* ------------------------------------------------------------------- */
190: /*
191:    FormInitialGuess - Computes initial guess.

193:    Input/Output Parameter:
194: .  x - the solution vector
195: */
196: PetscErrorCode FormInitialGuess(Vec x)
197: {
199:    PetscScalar    pfive = .50;
200:    VecSet(x,pfive);
201:    return 0;
202: }
203: /* ------------------------------------------------------------------- */
206: /* 
207:    FormFunction - Evaluates nonlinear function, F(x).

209:    Input Parameters:
210: .  snes - the SNES context
211: .  x - input vector
212: .  ctx - optional user-defined context, as set by SNESSetFunction()

214:    Output Parameter:
215: .  f - function vector

217:    Note:
218:    The user-defined context can contain any application-specific data
219:    needed for the function evaluation (such as various parameters, work
220:    vectors, and grid information).  In this program the context is just
221:    a vector containing the right-hand-side of the discretized PDE.
222:  */

224: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
225: {
226:    Vec            g = (Vec)ctx;
227:    PetscScalar    *xx,*ff,*gg,d;
229:    PetscInt       i,n;

231:   /*
232:      Get pointers to vector data.
233:        - For default PETSc vectors, VecGetArray() returns a pointer to
234:          the data array.  Otherwise, the routine is implementation dependent.
235:        - You MUST call VecRestoreArray() when you no longer need access to
236:          the array.
237:   */
238:    VecGetArray(x,&xx);
239:    VecGetArray(f,&ff);
240:    VecGetArray(g,&gg);

242:   /*
243:      Compute function
244:   */
245:    VecGetSize(x,&n);
246:    d = (PetscReal)(n - 1); d = d*d;
247:    ff[0]   = xx[0];
248:    for (i=1; i<n-1; i++) {
249:      ff[i] = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - gg[i];
250:    }
251:    ff[n-1] = xx[n-1] - 1.0;

253:   /*
254:      Restore vectors
255:   */
256:   VecRestoreArray(x,&xx);
257:   VecRestoreArray(f,&ff);
258:   VecRestoreArray(g,&gg);
259:   return 0;
260: }
261: /* ------------------------------------------------------------------- */
264: /*
265:    FormJacobian - Evaluates Jacobian matrix.

267:    Input Parameters:
268: .  snes - the SNES context
269: .  x - input vector
270: .  dummy - optional user-defined context (not used here)

272:    Output Parameters:
273: .  jac - Jacobian matrix
274: .  B - optionally different preconditioning matrix
275: .  flag - flag indicating matrix structure
276: */

278: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure*flag,void *dummy)
279: {
280:   PetscScalar    *xx,A[3],d;
282:   PetscInt       i,n,j[3];

284:   /*
285:      Get pointer to vector data
286:   */
287:   VecGetArray(x,&xx);

289:   /*
290:      Compute Jacobian entries and insert into matrix.
291:       - Note that in this case we set all elements for a particular
292:         row at once.
293:   */
294:   VecGetSize(x,&n);
295:   d = (PetscReal)(n - 1); d = d*d;

297:   /*
298:      Interior grid points
299:   */
300:   for (i=1; i<n-1; i++) {
301:     j[0] = i - 1; j[1] = i; j[2] = i + 1;
302:     A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
303:     MatSetValues(*jac,1,&i,3,j,A,INSERT_VALUES);
304:   }

306:   /*
307:      Boundary points
308:   */
309:   i = 0;   A[0] = 1.0;
310:   MatSetValues(*jac,1,&i,1,&i,A,INSERT_VALUES);
311:   i = n-1; A[0] = 1.0;
312:   MatSetValues(*jac,1,&i,1,&i,A,INSERT_VALUES);

314:   /*
315:      Restore vector
316:   */
317:   VecRestoreArray(x,&xx);

319:   /*
320:      Assemble matrix
321:   */
322:   MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
323:   MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);

325:   return 0;
326: }
327: /* ------------------------------------------------------------------- */
330: /*
331:    Monitor - User-defined monitoring routine that views the
332:    current iterate with an x-window plot.

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

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

351:   PetscPrintf(PETSC_COMM_WORLD,"iter = %D, SNES Function norm %G\n",its,fnorm);
352:   SNESGetSolution(snes,&x);
353:   VecView(x,monP->viewer);
354:   return 0;
355: }