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