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*/