Actual source code: ex6.c

petsc-main 2021-04-20
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 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 */
 37:   PetscInt       its,n = 5,i;
 38:   PetscMPIInt    size;
 39:   PetscScalar    h,xp,v;
 40:   MPI_Comm       comm;

 42:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 43:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 44:   if (size != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This is a uniprocessor example only!");
 45:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 46:   h    = 1.0/(n-1);
 47:   comm = PETSC_COMM_WORLD;
 48:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 49:      Create nonlinear solver context
 50:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 52:   SNESCreate(comm,&snes);

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

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


 73:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 74:      Create matrix data structure; set Jacobian evaluation routine
 75:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

 95:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 96:      Customize nonlinear solver; set runtime options
 97:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

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

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

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

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

183:    Output Parameter:
184: .  f - function vector

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

193: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
194: {
195:   Vec               g = (Vec)ctx;
196:   const PetscScalar *xx,*gg;
197:   PetscScalar       *ff,d;
198:   PetscErrorCode    ierr;
199:   PetscInt          i,n;

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

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

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

233:    Input Parameters:
234: .  snes - the SNES context
235: .  x - input vector
236: .  dummy - optional user-defined context (not used here)

238:    Output Parameters:
239: .  jac - Jacobian matrix
240: .  B - optionally different preconditioning matrix

242: */

244: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
245: {
246:   const PetscScalar *xx;
247:   PetscScalar       A[3],d;
248:   PetscErrorCode    ierr;
249:   PetscInt          i,n,j[3];

251:   /*
252:      Get pointer to vector data
253:   */
254:   VecGetArrayRead(x,&xx);

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

264:   /*
265:      Interior grid points
266:   */
267:   for (i=1; i<n-1; i++) {
268:     j[0] = i - 1; j[1] = i; j[2] = i + 1;
269:     A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
270:     MatSetValues(B,1,&i,3,j,A,INSERT_VALUES);
271:   }

273:   /*
274:      Boundary points
275:   */
276:   i = 0;   A[0] = 1.0;

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

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

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

284:   /*
285:      Restore vector
286:   */
287:   VecRestoreArrayRead(x,&xx);

289:   /*
290:      Assemble matrix
291:   */
292:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
293:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
294:   if (jac != B) {
295:     MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
296:     MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
297:   }
298:   return 0;
299: }

301: PetscErrorCode MySNESConvergedReasonView(SNES snes,void *ctx)
302: {
303:   PetscErrorCode        ierr;
304:   ReasonViewCtx         *monP = (ReasonViewCtx*) ctx;
305:   PetscViewer           viewer = monP->viewer;
306:   SNESConvergedReason   reason;
307:   const char            *strreason;

309:   SNESGetConvergedReason(snes,&reason);
310:   SNESGetConvergedReasonString(snes,&strreason);
311:   PetscViewerASCIIPrintf(viewer,"Customized SNES converged reason view\n");
312:   PetscViewerASCIIAddTab(viewer,1);
313:   if (reason > 0) {
314:     PetscViewerASCIIPrintf(viewer,"Converged due to %s\n",strreason);
315:   } else if (reason <= 0) {
316:     PetscViewerASCIIPrintf(viewer,"Did not converge due to %s\n",strreason);
317:   }
318:   PetscViewerASCIISubtractTab(viewer,1);
319:   return 0;
320: }

322: PetscErrorCode MyKSPConvergedReasonView(KSP ksp,void *ctx)
323: {
324:   PetscErrorCode        ierr;
325:   ReasonViewCtx         *monP = (ReasonViewCtx*) ctx;
326:   PetscViewer           viewer = monP->viewer;
327:   KSPConvergedReason    reason;
328:   const char            *reasonstr;

330:   KSPGetConvergedReason(ksp,&reason);
331:   KSPGetConvergedReasonString(ksp,&reasonstr);
332:   PetscViewerASCIIAddTab(viewer,2);
333:   PetscViewerASCIIPrintf(viewer,"Customized KSP converged reason view\n");
334:   PetscViewerASCIIAddTab(viewer,1);
335:   if (reason > 0) {
336:     PetscViewerASCIIPrintf(viewer,"Converged due to %s\n",reasonstr);
337:   } else if (reason <= 0) {
338:     PetscViewerASCIIPrintf(viewer,"Did not converge due to %s\n",reasonstr);
339:   }
340:   PetscViewerASCIISubtractTab(viewer,3);
341:   return 0;
342: }


345: /*TEST

347:    test:
348:       suffix: 1
349:       nsize: 1

351:    test:
352:       suffix: 2
353:       nsize: 1
354:       args: -ksp_converged_reason_view_cancel

356:    test:
357:       suffix: 3
358:       nsize: 1
359:       args: -ksp_converged_reason_view_cancel -ksp_converged_reason

361:    test:
362:       suffix: 4
363:       nsize: 1
364:       args: -snes_converged_reason_view_cancel

366:    test:
367:       suffix: 5
368:       nsize: 1
369:       args: -snes_converged_reason_view_cancel -snes_converged_reason

371: TEST*/