Actual source code: ex6.c


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

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

 11: #include <petscsnes.h>

 13: /*
 14:    User-defined routines
 15: */
 16: PETSC_EXTERN PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
 17: PETSC_EXTERN PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
 18: PETSC_EXTERN PetscErrorCode FormInitialGuess(Vec);
 19: PETSC_EXTERN PetscErrorCode MySNESConvergedReasonView(SNES,void*);
 20: PETSC_EXTERN PetscErrorCode MyKSPConvergedReasonView(KSP,void*);

 22: /*
 23:    User-defined context for monitoring
 24: */
 25: typedef struct {
 26:   PetscViewer viewer;
 27: } ReasonViewCtx;

 29: int main(int argc,char **argv)
 30: {
 31:   SNES           snes;                /* SNES context */
 32:   KSP            ksp;                 /* KSP context */
 33:   Vec            x,r,F,U;             /* vectors */
 34:   Mat            J;                   /* Jacobian matrix */
 35:   ReasonViewCtx     monP;             /* monitoring context */
 36:   PetscInt       its,n = 5,i;
 37:   PetscMPIInt    size;
 38:   PetscScalar    h,xp,v;
 39:   MPI_Comm       comm;

 41:   PetscInitialize(&argc,&argv,(char*)0,help);
 42:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 44:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 45:   h    = 1.0/(n-1);
 46:   comm = PETSC_COMM_WORLD;
 47:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 48:      Create nonlinear solver context
 49:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 51:   SNESCreate(comm,&snes);

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

 66:   /*
 67:      Set function evaluation routine and vector
 68:   */
 69:   SNESSetFunction(snes,r,FormFunction,(void*)F);

 71:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 72:      Create matrix data structure; set Jacobian evaluation routine
 73:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 75:   MatCreate(comm,&J);
 76:   MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);
 77:   MatSetFromOptions(J);
 78:   MatSeqAIJSetPreallocation(J,3,NULL);

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

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

 93:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 94:      Customize nonlinear solver; set runtime options
 95:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 97:   /*
 98:      Set an optional user-defined reasonview routine
 99:   */
100:   PetscViewerASCIIGetStdout(comm,&monP.viewer);
101:   /* Just make sure we can not repeat addding the same function
102:    * PETSc will be able to igore the repeated function
103:    */
104:   for (i=0; i<4; i++) {
105:     SNESConvergedReasonViewSet(snes,MySNESConvergedReasonView,&monP,0);
106:   }
107:   SNESGetKSP(snes,&ksp);
108:   /* Just make sure we can not repeat addding the same function
109:    * PETSc will be able to igore the repeated function
110:    */
111:   for (i=0; i<4; i++) {
112:     KSPConvergedReasonViewSet(ksp,MyKSPConvergedReasonView,&monP,0);
113:   }
114:   /*
115:      Set SNES/KSP/KSP/PC runtime options, e.g.,
116:          -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
117:   */
118:   SNESSetFromOptions(snes);

120:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121:      Initialize application:
122:      Store right-hand-side of PDE and exact solution
123:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

125:   xp = 0.0;
126:   for (i=0; i<n; i++) {
127:     v    = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
128:     VecSetValues(F,1,&i,&v,INSERT_VALUES);
129:     v    = xp*xp*xp;
130:     VecSetValues(U,1,&i,&v,INSERT_VALUES);
131:     xp  += h;
132:   }

134:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135:      Evaluate initial guess; then solve nonlinear system
136:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
137:   /*
138:      Note: The user should initialize the vector, x, with the initial guess
139:      for the nonlinear solver prior to calling SNESSolve().  In particular,
140:      to employ an initial guess of zero, the user should explicitly set
141:      this vector to zero by calling VecSet().
142:   */
143:   FormInitialGuess(x);
144:   SNESSolve(snes,NULL,x);
145:   SNESGetIterationNumber(snes,&its);

147:   /*
148:      Free work space.  All PETSc objects should be destroyed when they
149:      are no longer needed.
150:   */
151:   VecDestroy(&x));  PetscCall(VecDestroy(&r);
152:   VecDestroy(&U));  PetscCall(VecDestroy(&F);
153:   MatDestroy(&J));  PetscCall(SNESDestroy(&snes);
154:   /*PetscViewerDestroy(&monP.viewer);*/
155:   PetscFinalize();
156:   return 0;
157: }
158: /* ------------------------------------------------------------------- */
159: /*
160:    FormInitialGuess - Computes initial guess.

162:    Input/Output Parameter:
163: .  x - the solution vector
164: */
165: PetscErrorCode FormInitialGuess(Vec x)
166: {
167:   PetscScalar    pfive = .50;
168:   VecSet(x,pfive);
169:   return 0;
170: }
171: /* ------------------------------------------------------------------- */
172: /*
173:    FormFunction - Evaluates nonlinear function, F(x).

175:    Input Parameters:
176: .  snes - the SNES context
177: .  x - input vector
178: .  ctx - optional user-defined context, as set by SNESSetFunction()

180:    Output Parameter:
181: .  f - function vector

183:    Note:
184:    The user-defined context can contain any application-specific data
185:    needed for the function evaluation (such as various parameters, work
186:    vectors, and grid information).  In this program the context is just
187:    a vector containing the right-hand-side of the discretized PDE.
188:  */

190: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
191: {
192:   Vec               g = (Vec)ctx;
193:   const PetscScalar *xx,*gg;
194:   PetscScalar       *ff,d;
195:   PetscInt          i,n;

197:   /*
198:      Get pointers to vector data.
199:        - For default PETSc vectors, VecGetArray() returns a pointer to
200:          the data array.  Otherwise, the routine is implementation dependent.
201:        - You MUST call VecRestoreArray() when you no longer need access to
202:          the array.
203:   */
204:   VecGetArrayRead(x,&xx);
205:   VecGetArray(f,&ff);
206:   VecGetArrayRead(g,&gg);

208:   /*
209:      Compute function
210:   */
211:   VecGetSize(x,&n);
212:   d     = (PetscReal)(n - 1); d = d*d;
213:   ff[0] = xx[0];
214:   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];
215:   ff[n-1] = xx[n-1] - 1.0;

217:   /*
218:      Restore vectors
219:   */
220:   VecRestoreArrayRead(x,&xx);
221:   VecRestoreArray(f,&ff);
222:   VecRestoreArrayRead(g,&gg);
223:   return 0;
224: }
225: /* ------------------------------------------------------------------- */
226: /*
227:    FormJacobian - Evaluates Jacobian matrix.

229:    Input Parameters:
230: .  snes - the SNES context
231: .  x - input vector
232: .  dummy - optional user-defined context (not used here)

234:    Output Parameters:
235: .  jac - Jacobian matrix
236: .  B - optionally different preconditioning matrix

238: */

240: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
241: {
242:   const PetscScalar *xx;
243:   PetscScalar       A[3],d;
244:   PetscInt          i,n,j[3];

246:   /*
247:      Get pointer to vector data
248:   */
249:   VecGetArrayRead(x,&xx);

251:   /*
252:      Compute Jacobian entries and insert into matrix.
253:       - Note that in this case we set all elements for a particular
254:         row at once.
255:   */
256:   VecGetSize(x,&n);
257:   d    = (PetscReal)(n - 1); d = d*d;

259:   /*
260:      Interior grid points
261:   */
262:   for (i=1; i<n-1; i++) {
263:     j[0] = i - 1; j[1] = i; j[2] = i + 1;
264:     A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
265:     MatSetValues(B,1,&i,3,j,A,INSERT_VALUES);
266:   }

268:   /*
269:      Boundary points
270:   */
271:   i = 0;   A[0] = 1.0;

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

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

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

279:   /*
280:      Restore vector
281:   */
282:   VecRestoreArrayRead(x,&xx);

284:   /*
285:      Assemble matrix
286:   */
287:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
288:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
289:   if (jac != B) {
290:     MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
291:     MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
292:   }
293:   return 0;
294: }

296: PetscErrorCode MySNESConvergedReasonView(SNES snes,void *ctx)
297: {
298:   ReasonViewCtx         *monP = (ReasonViewCtx*) ctx;
299:   PetscViewer           viewer = monP->viewer;
300:   SNESConvergedReason   reason;
301:   const char            *strreason;

303:   SNESGetConvergedReason(snes,&reason);
304:   SNESGetConvergedReasonString(snes,&strreason);
305:   PetscViewerASCIIPrintf(viewer,"Customized SNES converged reason view\n");
306:   PetscViewerASCIIAddTab(viewer,1);
307:   if (reason > 0) {
308:     PetscViewerASCIIPrintf(viewer,"Converged due to %s\n",strreason);
309:   } else if (reason <= 0) {
310:     PetscViewerASCIIPrintf(viewer,"Did not converge due to %s\n",strreason);
311:   }
312:   PetscViewerASCIISubtractTab(viewer,1);
313:   return 0;
314: }

316: PetscErrorCode MyKSPConvergedReasonView(KSP ksp,void *ctx)
317: {
318:   ReasonViewCtx         *monP = (ReasonViewCtx*) ctx;
319:   PetscViewer           viewer = monP->viewer;
320:   KSPConvergedReason    reason;
321:   const char            *reasonstr;

323:   KSPGetConvergedReason(ksp,&reason);
324:   KSPGetConvergedReasonString(ksp,&reasonstr);
325:   PetscViewerASCIIAddTab(viewer,2);
326:   PetscViewerASCIIPrintf(viewer,"Customized KSP converged reason view\n");
327:   PetscViewerASCIIAddTab(viewer,1);
328:   if (reason > 0) {
329:     PetscViewerASCIIPrintf(viewer,"Converged due to %s\n",reasonstr);
330:   } else if (reason <= 0) {
331:     PetscViewerASCIIPrintf(viewer,"Did not converge due to %s\n",reasonstr);
332:   }
333:   PetscViewerASCIISubtractTab(viewer,3);
334:   return 0;
335: }

337: /*TEST

339:    test:
340:       suffix: 1
341:       nsize: 1
342:       filter: sed -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g"

344:    test:
345:       suffix: 2
346:       nsize: 1
347:       args: -ksp_converged_reason_view_cancel
348:       filter: sed -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g"

350:    test:
351:       suffix: 3
352:       nsize: 1
353:       args: -ksp_converged_reason_view_cancel -ksp_converged_reason
354:       filter: sed -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g"

356:    test:
357:       suffix: 4
358:       nsize: 1
359:       args: -snes_converged_reason_view_cancel
360:       filter: sed -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g"

362:    test:
363:       suffix: 5
364:       nsize: 1
365:       args: -snes_converged_reason_view_cancel -snes_converged_reason
366:       filter: sed -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g"

368: TEST*/