Actual source code: ex2.c
2: static char help[] = "Newton method to solve u'' + u^{2} = f, sequentially.\n\
3: This example employs a user-defined monitoring routine.\n\n";
5: /*T
6: Concepts: SNES^basic uniprocessor example
7: Concepts: SNES^setting a user-defined monitoring routine
8: Processors: 1
9: T*/
11: /*
12: Include "petscdraw.h" so that we can use PETSc drawing routines.
13: Include "petscsnes.h" so that we can use SNES solvers. Note that this
14: file automatically includes:
15: petscsys.h - base PETSc routines petscvec.h - vectors
16: petscmat.h - matrices
17: petscis.h - index sets petscksp.h - Krylov subspace methods
18: petscviewer.h - viewers petscpc.h - preconditioners
19: petscksp.h - linear solvers
20: */
22: #include <petscsnes.h>
24: /*
25: User-defined routines
26: */
27: extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
28: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
29: extern PetscErrorCode FormInitialGuess(Vec);
30: extern PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*);
32: /*
33: User-defined context for monitoring
34: */
35: typedef struct {
36: PetscViewer viewer;
37: } MonitorCtx;
39: int main(int argc,char **argv)
40: {
41: SNES snes; /* SNES context */
42: Vec x,r,F,U; /* vectors */
43: Mat J; /* Jacobian matrix */
44: MonitorCtx monP; /* monitoring context */
45: PetscInt its,n = 5,i,maxit,maxf;
46: PetscMPIInt size;
47: PetscScalar h,xp,v,none = -1.0;
48: PetscReal abstol,rtol,stol,norm;
50: PetscInitialize(&argc,&argv,(char*)0,help);
51: MPI_Comm_size(PETSC_COMM_WORLD,&size);
53: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
54: h = 1.0/(n-1);
56: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
57: Create nonlinear solver context
58: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
60: SNESCreate(PETSC_COMM_WORLD,&snes);
62: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
63: Create vector data structures; set function evaluation routine
64: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
65: /*
66: Note that we form 1 vector from scratch and then duplicate as needed.
67: */
68: VecCreate(PETSC_COMM_WORLD,&x);
69: VecSetSizes(x,PETSC_DECIDE,n);
70: VecSetFromOptions(x);
71: VecDuplicate(x,&r);
72: VecDuplicate(x,&F);
73: VecDuplicate(x,&U);
75: /*
76: Set function evaluation routine and vector
77: */
78: SNESSetFunction(snes,r,FormFunction,(void*)F);
80: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
81: Create matrix data structure; set Jacobian evaluation routine
82: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
84: MatCreate(PETSC_COMM_WORLD,&J);
85: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);
86: MatSetFromOptions(J);
87: MatSeqAIJSetPreallocation(J,3,NULL);
89: /*
90: Set Jacobian matrix data structure and default Jacobian evaluation
91: routine. User can override with:
92: -snes_fd : default finite differencing approximation of Jacobian
93: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
94: (unless user explicitly sets preconditioner)
95: -snes_mf_operator : form preconditioning matrix as set by the user,
96: but use matrix-free approx for Jacobian-vector
97: products within Newton-Krylov method
98: */
100: SNESSetJacobian(snes,J,J,FormJacobian,NULL);
102: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103: Customize nonlinear solver; set runtime options
104: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
106: /*
107: Set an optional user-defined monitoring routine
108: */
109: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,400,400,&monP.viewer);
110: SNESMonitorSet(snes,Monitor,&monP,0);
112: /*
113: Set names for some vectors to facilitate monitoring (optional)
114: */
115: PetscObjectSetName((PetscObject)x,"Approximate Solution");
116: PetscObjectSetName((PetscObject)U,"Exact Solution");
118: /*
119: Set SNES/KSP/KSP/PC runtime options, e.g.,
120: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
121: */
122: SNESSetFromOptions(snes);
124: /*
125: Print parameters used for convergence testing (optional) ... just
126: to demonstrate this routine; this information is also printed with
127: the option -snes_view
128: */
129: SNESGetTolerances(snes,&abstol,&rtol,&stol,&maxit,&maxf);
130: PetscPrintf(PETSC_COMM_WORLD,"atol=%g, rtol=%g, stol=%g, maxit=%D, maxf=%D\n",(double)abstol,(double)rtol,(double)stol,maxit,maxf);
132: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
133: Initialize application:
134: Store right-hand-side of PDE and exact solution
135: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
137: xp = 0.0;
138: for (i=0; i<n; i++) {
139: v = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
140: VecSetValues(F,1,&i,&v,INSERT_VALUES);
141: v = xp*xp*xp;
142: VecSetValues(U,1,&i,&v,INSERT_VALUES);
143: xp += h;
144: }
146: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: Evaluate initial guess; then solve nonlinear system
148: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
149: /*
150: Note: The user should initialize the vector, x, with the initial guess
151: for the nonlinear solver prior to calling SNESSolve(). In particular,
152: to employ an initial guess of zero, the user should explicitly set
153: this vector to zero by calling VecSet().
154: */
155: FormInitialGuess(x);
156: SNESSolve(snes,NULL,x);
157: SNESGetIterationNumber(snes,&its);
158: PetscPrintf(PETSC_COMM_WORLD,"number of SNES iterations = %D\n\n",its);
160: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
161: Check solution and clean up
162: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
164: /*
165: Check the error
166: */
167: VecAXPY(x,none,U);
168: VecNorm(x,NORM_2,&norm);
169: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);
171: /*
172: Free work space. All PETSc objects should be destroyed when they
173: are no longer needed.
174: */
175: VecDestroy(&x)); PetscCall(VecDestroy(&r);
176: VecDestroy(&U)); PetscCall(VecDestroy(&F);
177: MatDestroy(&J)); PetscCall(SNESDestroy(&snes);
178: PetscViewerDestroy(&monP.viewer);
179: PetscFinalize();
180: return 0;
181: }
182: /* ------------------------------------------------------------------- */
183: /*
184: FormInitialGuess - Computes initial guess.
186: Input/Output Parameter:
187: . x - the solution vector
188: */
189: PetscErrorCode FormInitialGuess(Vec x)
190: {
191: PetscScalar pfive = .50;
192: VecSet(x,pfive);
193: return 0;
194: }
195: /* ------------------------------------------------------------------- */
196: /*
197: FormFunction - Evaluates nonlinear function, F(x).
199: Input Parameters:
200: . snes - the SNES context
201: . x - input vector
202: . ctx - optional user-defined context, as set by SNESSetFunction()
204: Output Parameter:
205: . f - function vector
207: Note:
208: The user-defined context can contain any application-specific data
209: needed for the function evaluation (such as various parameters, work
210: vectors, and grid information). In this program the context is just
211: a vector containing the right-hand-side of the discretized PDE.
212: */
214: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
215: {
216: Vec g = (Vec)ctx;
217: const PetscScalar *xx,*gg;
218: PetscScalar *ff,d;
219: PetscInt i,n;
221: /*
222: Get pointers to vector data.
223: - For default PETSc vectors, VecGetArray() returns a pointer to
224: the data array. Otherwise, the routine is implementation dependent.
225: - You MUST call VecRestoreArray() when you no longer need access to
226: the array.
227: */
228: VecGetArrayRead(x,&xx);
229: VecGetArray(f,&ff);
230: VecGetArrayRead(g,&gg);
232: /*
233: Compute function
234: */
235: VecGetSize(x,&n);
236: d = (PetscReal)(n - 1); d = d*d;
237: ff[0] = xx[0];
238: 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];
239: ff[n-1] = xx[n-1] - 1.0;
241: /*
242: Restore vectors
243: */
244: VecRestoreArrayRead(x,&xx);
245: VecRestoreArray(f,&ff);
246: VecRestoreArrayRead(g,&gg);
247: return 0;
248: }
249: /* ------------------------------------------------------------------- */
250: /*
251: FormJacobian - Evaluates Jacobian matrix.
253: Input Parameters:
254: . snes - the SNES context
255: . x - input vector
256: . dummy - optional user-defined context (not used here)
258: Output Parameters:
259: . jac - Jacobian matrix
260: . B - optionally different preconditioning matrix
262: */
264: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
265: {
266: const PetscScalar *xx;
267: PetscScalar A[3],d;
268: PetscInt i,n,j[3];
270: /*
271: Get pointer to vector data
272: */
273: VecGetArrayRead(x,&xx);
275: /*
276: Compute Jacobian entries and insert into matrix.
277: - Note that in this case we set all elements for a particular
278: row at once.
279: */
280: VecGetSize(x,&n);
281: d = (PetscReal)(n - 1); d = d*d;
283: /*
284: Interior grid points
285: */
286: for (i=1; i<n-1; i++) {
287: j[0] = i - 1; j[1] = i; j[2] = i + 1;
288: A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
289: MatSetValues(B,1,&i,3,j,A,INSERT_VALUES);
290: }
292: /*
293: Boundary points
294: */
295: i = 0; A[0] = 1.0;
297: MatSetValues(B,1,&i,1,&i,A,INSERT_VALUES);
299: i = n-1; A[0] = 1.0;
301: MatSetValues(B,1,&i,1,&i,A,INSERT_VALUES);
303: /*
304: Restore vector
305: */
306: VecRestoreArrayRead(x,&xx);
308: /*
309: Assemble matrix
310: */
311: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
312: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
313: if (jac != B) {
314: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
315: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
316: }
317: return 0;
318: }
319: /* ------------------------------------------------------------------- */
320: /*
321: Monitor - User-defined monitoring routine that views the
322: current iterate with an x-window plot.
324: Input Parameters:
325: snes - the SNES context
326: its - iteration number
327: norm - 2-norm function value (may be estimated)
328: ctx - optional user-defined context for private data for the
329: monitor routine, as set by SNESMonitorSet()
331: Note:
332: See the manpage for PetscViewerDrawOpen() for useful runtime options,
333: such as -nox to deactivate all x-window output.
334: */
335: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
336: {
337: MonitorCtx *monP = (MonitorCtx*) ctx;
338: Vec x;
340: PetscPrintf(PETSC_COMM_WORLD,"iter = %D, SNES Function norm %g\n",its,(double)fnorm);
341: SNESGetSolution(snes,&x);
342: VecView(x,monP->viewer);
343: return 0;
344: }
346: /*TEST
348: test:
349: args: -nox -snes_monitor_cancel -snes_monitor_short -snes_view -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always
351: test:
352: suffix: 2
353: args: -nox -snes_monitor_cancel -snes_monitor_short -snes_type newtontr -snes_view
354: requires: !single
356: test:
357: suffix: 3
358: args: -nox -malloc no -options_left no -snes_monitor_cancel -snes_monitor_short -snes_view -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always
360: test:
361: suffix: 4
362: args: -nox -snes_monitor_cancel -snes_monitor_short -snes_type newtontrdc -snes_view
363: requires: !single
365: TEST*/