Actual source code: ex2.c

petsc-3.9.3 2018-07-02
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*/



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

 24:  #include <petscsnes.h>

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

 34: /*
 35:    User-defined context for monitoring
 36: */
 37: typedef struct {
 38:   PetscViewer viewer;
 39: } 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);if (ierr) return ierr;
 54:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 55:   if (size != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This is a uniprocessor example only!");
 56:   PetscOptionsGetInt(NULL,NULL,"-n",&n,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);
 91:   MatSeqAIJSetPreallocation(J,3,NULL);

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

210:    Output Parameter:
211: .  f - function vector

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

220: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
221: {
222:   Vec               g = (Vec)ctx;
223:   const PetscScalar *xx,*gg;
224:   PetscScalar       *ff,d;
225:   PetscErrorCode    ierr;
226:   PetscInt          i,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:   VecGetArrayRead(x,&xx);
236:   VecGetArray(f,&ff);
237:   VecGetArrayRead(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++) ff[i] = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - gg[i];
246:   ff[n-1] = xx[n-1] - 1.0;

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

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

265:    Output Parameters:
266: .  jac - Jacobian matrix
267: .  B - optionally different preconditioning matrix

269: */

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

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

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

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

300:   /*
301:      Boundary points
302:   */
303:   i = 0;   A[0] = 1.0;

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

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

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

311:   /*
312:      Restore vector
313:   */
314:   VecRestoreArrayRead(x,&xx);

316:   /*
317:      Assemble matrix
318:   */
319:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
320:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
321:   if (jac != B) {
322:     MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
323:     MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
324:   }
325:   return 0;
326: }
327: /* ------------------------------------------------------------------- */
328: /*
329:    Monitor - User-defined monitoring routine that views the
330:    current iterate with an x-window plot.

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

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

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


356: /*TEST

358:    test:
359:       args: -nox -snes_monitor_cancel -snes_monitor_short -snes_view -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always

361:    test:
362:       suffix: 2
363:       args: -nox -snes_monitor_cancel -snes_monitor_short -snes_type newtontr -snes_view
364:       requires: !single

366:    test:
367:       suffix: 3
368:       args: -nox -snes_monitor_cancel -snes_monitor_short -malloc no -snes_view -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always

370: TEST*/