Actual source code: ex2.c

petsc-master 2017-06-22
Report Typos and Errors

  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 */
 46:   PetscInt       its,n = 5,i,maxit,maxf;
 47:   PetscMPIInt    size;
 48:   PetscScalar    h,xp,v,none = -1.0;
 49:   PetscReal      abstol,rtol,stol,norm;

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

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

 61:   SNESCreate(PETSC_COMM_WORLD,&snes);

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

203:    Input Parameters:
204: .  snes - the SNES context
205: .  x - input vector
206: .  ctx - optional user-defined context, as set by SNESSetFunction()

208:    Output Parameter:
209: .  f - function vector

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

218: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
219: {
220:   Vec               g = (Vec)ctx;
221:   const PetscScalar *xx,*gg;
222:   PetscScalar       *ff,d;
223:   PetscErrorCode    ierr;
224:   PetscInt          i,n;

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

237:   /*
238:      Compute function
239:   */
240:   VecGetSize(x,&n);
241:   d     = (PetscReal)(n - 1); d = d*d;
242:   ff[0] = xx[0];
243:   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];
244:   ff[n-1] = xx[n-1] - 1.0;

246:   /*
247:      Restore vectors
248:   */
249:   VecRestoreArrayRead(x,&xx);
250:   VecRestoreArray(f,&ff);
251:   VecRestoreArrayRead(g,&gg);
252:   return 0;
253: }
254: /* ------------------------------------------------------------------- */
255: /*
256:    FormJacobian - Evaluates Jacobian matrix.

258:    Input Parameters:
259: .  snes - the SNES context
260: .  x - input vector
261: .  dummy - optional user-defined context (not used here)

263:    Output Parameters:
264: .  jac - Jacobian matrix
265: .  B - optionally different preconditioning matrix

267: */

269: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
270: {
271:   const PetscScalar *xx;
272:   PetscScalar       A[3],d;
273:   PetscErrorCode    ierr;
274:   PetscInt          i,n,j[3];

276:   /*
277:      Get pointer to vector data
278:   */
279:   VecGetArrayRead(x,&xx);

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

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

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

303:   MatSetValues(jac,1,&i,1,&i,A,INSERT_VALUES);

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

307:   MatSetValues(jac,1,&i,1,&i,A,INSERT_VALUES);

309:   /*
310:      Restore vector
311:   */
312:   VecRestoreArrayRead(x,&xx);

314:   /*
315:      Assemble matrix
316:   */
317:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
318:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
319:   if (jac != B) {
320:     MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
321:     MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
322:   }
323:   return 0;
324: }
325: /* ------------------------------------------------------------------- */
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 SNESMonitorSet()

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

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