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: #include <petscsnes.h>

  7: /*
  8:    User-defined routines
  9: */
 10: PETSC_EXTERN PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
 11: PETSC_EXTERN PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
 12: PETSC_EXTERN PetscErrorCode FormInitialGuess(Vec);
 13: PETSC_EXTERN PetscErrorCode MySNESConvergedReasonView(SNES,void*);
 14: PETSC_EXTERN PetscErrorCode MyKSPConvergedReasonView(KSP,void*);

 16: /*
 17:    User-defined context for monitoring
 18: */
 19: typedef struct {
 20:   PetscViewer viewer;
 21: } ReasonViewCtx;

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

 35:   PetscInitialize(&argc,&argv,(char*)0,help);
 36:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 38:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 39:   h    = 1.0/(n-1);
 40:   comm = PETSC_COMM_WORLD;
 41:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 42:      Create nonlinear solver context
 43:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 45:   SNESCreate(comm,&snes);

 47:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 48:      Create vector data structures; set function evaluation routine
 49:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 50:   /*
 51:      Note that we form 1 vector from scratch and then duplicate as needed.
 52:   */
 53:   VecCreate(comm,&x);
 54:   VecSetSizes(x,PETSC_DECIDE,n);
 55:   VecSetFromOptions(x);
 56:   VecDuplicate(x,&r);
 57:   VecDuplicate(x,&F);
 58:   VecDuplicate(x,&U);

 60:   /*
 61:      Set function evaluation routine and vector
 62:   */
 63:   SNESSetFunction(snes,r,FormFunction,(void*)F);

 65:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 66:      Create matrix data structure; set Jacobian evaluation routine
 67:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 69:   MatCreate(comm,&J);
 70:   MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);
 71:   MatSetFromOptions(J);
 72:   MatSeqAIJSetPreallocation(J,3,NULL);

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

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

 87:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 88:      Customize nonlinear solver; set runtime options
 89:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

114:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115:      Initialize application:
116:      Store right-hand-side of PDE and exact solution
117:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

119:   xp = 0.0;
120:   for (i=0; i<n; i++) {
121:     v    = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
122:     VecSetValues(F,1,&i,&v,INSERT_VALUES);
123:     v    = xp*xp*xp;
124:     VecSetValues(U,1,&i,&v,INSERT_VALUES);
125:     xp  += h;
126:   }

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

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

156:    Input/Output Parameter:
157: .  x - the solution vector
158: */
159: PetscErrorCode FormInitialGuess(Vec x)
160: {
161:   PetscScalar    pfive = .50;
162:   VecSet(x,pfive);
163:   return 0;
164: }
165: /* ------------------------------------------------------------------- */
166: /*
167:    FormFunction - Evaluates nonlinear function, F(x).

169:    Input Parameters:
170: .  snes - the SNES context
171: .  x - input vector
172: .  ctx - optional user-defined context, as set by SNESSetFunction()

174:    Output Parameter:
175: .  f - function vector

177:    Note:
178:    The user-defined context can contain any application-specific data
179:    needed for the function evaluation (such as various parameters, work
180:    vectors, and grid information).  In this program the context is just
181:    a vector containing the right-hand-side of the discretized PDE.
182:  */

184: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
185: {
186:   Vec               g = (Vec)ctx;
187:   const PetscScalar *xx,*gg;
188:   PetscScalar       *ff,d;
189:   PetscInt          i,n;

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

202:   /*
203:      Compute function
204:   */
205:   VecGetSize(x,&n);
206:   d     = (PetscReal)(n - 1); d = d*d;
207:   ff[0] = xx[0];
208:   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];
209:   ff[n-1] = xx[n-1] - 1.0;

211:   /*
212:      Restore vectors
213:   */
214:   VecRestoreArrayRead(x,&xx);
215:   VecRestoreArray(f,&ff);
216:   VecRestoreArrayRead(g,&gg);
217:   return 0;
218: }
219: /* ------------------------------------------------------------------- */
220: /*
221:    FormJacobian - Evaluates Jacobian matrix.

223:    Input Parameters:
224: .  snes - the SNES context
225: .  x - input vector
226: .  dummy - optional user-defined context (not used here)

228:    Output Parameters:
229: .  jac - Jacobian matrix
230: .  B - optionally different preconditioning matrix

232: */

234: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
235: {
236:   const PetscScalar *xx;
237:   PetscScalar       A[3],d;
238:   PetscInt          i,n,j[3];

240:   /*
241:      Get pointer to vector data
242:   */
243:   VecGetArrayRead(x,&xx);

245:   /*
246:      Compute Jacobian entries and insert into matrix.
247:       - Note that in this case we set all elements for a particular
248:         row at once.
249:   */
250:   VecGetSize(x,&n);
251:   d    = (PetscReal)(n - 1); d = d*d;

253:   /*
254:      Interior grid points
255:   */
256:   for (i=1; i<n-1; i++) {
257:     j[0] = i - 1; j[1] = i; j[2] = i + 1;
258:     A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
259:     MatSetValues(B,1,&i,3,j,A,INSERT_VALUES);
260:   }

262:   /*
263:      Boundary points
264:   */
265:   i = 0;   A[0] = 1.0;

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

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

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

273:   /*
274:      Restore vector
275:   */
276:   VecRestoreArrayRead(x,&xx);

278:   /*
279:      Assemble matrix
280:   */
281:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
282:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
283:   if (jac != B) {
284:     MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
285:     MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
286:   }
287:   return 0;
288: }

290: PetscErrorCode MySNESConvergedReasonView(SNES snes,void *ctx)
291: {
292:   ReasonViewCtx         *monP = (ReasonViewCtx*) ctx;
293:   PetscViewer           viewer = monP->viewer;
294:   SNESConvergedReason   reason;
295:   const char            *strreason;

297:   SNESGetConvergedReason(snes,&reason);
298:   SNESGetConvergedReasonString(snes,&strreason);
299:   PetscViewerASCIIPrintf(viewer,"Customized SNES converged reason view\n");
300:   PetscViewerASCIIAddTab(viewer,1);
301:   if (reason > 0) {
302:     PetscViewerASCIIPrintf(viewer,"Converged due to %s\n",strreason);
303:   } else if (reason <= 0) {
304:     PetscViewerASCIIPrintf(viewer,"Did not converge due to %s\n",strreason);
305:   }
306:   PetscViewerASCIISubtractTab(viewer,1);
307:   return 0;
308: }

310: PetscErrorCode MyKSPConvergedReasonView(KSP ksp,void *ctx)
311: {
312:   ReasonViewCtx         *monP = (ReasonViewCtx*) ctx;
313:   PetscViewer           viewer = monP->viewer;
314:   KSPConvergedReason    reason;
315:   const char            *reasonstr;

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

331: /*TEST

333:    test:
334:       suffix: 1
335:       nsize: 1
336:       filter: sed -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g"

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

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

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

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

362: TEST*/