Actual source code: ex2.c

  1: /*$Id: ex2.c,v 1.83 2001/08/07 21:31:17 bsmith Exp $*/

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

  6: /*T
  7:    Concepts: SNES^basic uniprocessor example
  8:    Concepts: SNES^setting a user-defined monitoring routine
  9:    Processors: 1
 10: T*/

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

 23:  #include petscsnes.h

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

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

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

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

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

 62:   SNESCreate(PETSC_COMM_WORLD,&snes);

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

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


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

 87:   MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,&J);
 88:   MatSetFromOptions(J);

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

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

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

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

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

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

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

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

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

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

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

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


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

183:   return 0;
184: }
185: /* ------------------------------------------------------------------- */
188: /*
189:    FormInitialGuess - Computes initial guess.

191:    Input/Output Parameter:
192: .  x - the solution vector
193: */
194: int FormInitialGuess(Vec x)
195: {
196:    int    ierr;
197:    PetscScalar pfive = .50;
198:    VecSet(&pfive,x);
199:    return 0;
200: }
201: /* ------------------------------------------------------------------- */
204: /* 
205:    FormFunction - Evaluates nonlinear function, F(x).

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

212:    Output Parameter:
213: .  f - function vector

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

222: int FormFunction(SNES snes,Vec x,Vec f,void *ctx)
223: {
224:    Vec    g = (Vec)ctx;
225:    PetscScalar *xx,*ff,*gg,d;
226:    int    i,ierr,n;

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

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

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

264:    Input Parameters:
265: .  snes - the SNES context
266: .  x - input vector
267: .  dummy - optional user-defined context (not used here)

269:    Output Parameters:
270: .  jac - Jacobian matrix
271: .  B - optionally different preconditioning matrix
272: .  flag - flag indicating matrix structure
273: */

275: int FormJacobian(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure*flag,void *dummy)
276: {
277:   PetscScalar *xx,A[3],d;
278:   int    i,n,j[3],ierr;

280:   /*
281:      Get pointer to vector data
282:   */
283:   VecGetArray(x,&xx);

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

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

302:   /*
303:      Boundary points
304:   */
305:   i = 0;   A[0] = 1.0;
306:   MatSetValues(*jac,1,&i,1,&i,A,INSERT_VALUES);
307:   i = n-1; A[0] = 1.0;
308:   MatSetValues(*jac,1,&i,1,&i,A,INSERT_VALUES);

310:   /*
311:      Restore vector
312:   */
313:   VecRestoreArray(x,&xx);

315:   /*
316:      Assemble matrix
317:   */
318:   MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
319:   MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);

321:   return 0;
322: }
323: /* ------------------------------------------------------------------- */
326: /*
327:    Monitor - User-defined monitoring routine that views the
328:    current iterate with an x-window plot.

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

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

347:   PetscPrintf(PETSC_COMM_WORLD,"iter = %d, SNES Function norm %g\n",its,fnorm);
348:   SNESGetSolution(snes,&x);
349:   VecView(x,monP->viewer);
350:   return 0;
351: }