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: petsc.h - base PETSc routines petscvec.h - vectors
16: petscsys.h - system routines 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: */
32: /*
33: User-defined context for monitoring
34: */
35: typedef struct {
36: PetscViewer viewer;
37: } MonitorCtx;
41: int main(int argc,char **argv)
42: {
43: SNES snes; /* SNES context */
44: Vec x,r,F,U; /* vectors */
45: Mat J; /* Jacobian matrix */
46: MonitorCtx monP; /* monitoring context */
48: PetscInt its,n = 5,i,maxit,maxf;
49: PetscMPIInt size;
50: PetscScalar h,xp,v,none = -1.0;
51: PetscReal abstol,rtol,stol,norm;
53: PetscInitialize(&argc,&argv,(char *)0,help);
54: MPI_Comm_size(PETSC_COMM_WORLD,&size);
55: if (size != 1) SETERRQ(1,"This is a uniprocessor example only!");
56: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
57: h = 1.0/(n-1);
59: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
60: Create nonlinear solver context
61: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
63: SNESCreate(PETSC_COMM_WORLD,&snes);
65: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
66: Create vector data structures; set function evaluation routine
67: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
68: /*
69: Note that we form 1 vector from scratch and then duplicate as needed.
70: */
71: VecCreate(PETSC_COMM_WORLD,&x);
72: VecSetSizes(x,PETSC_DECIDE,n);
73: VecSetFromOptions(x);
74: VecDuplicate(x,&r);
75: VecDuplicate(x,&F);
76: VecDuplicate(x,&U);
78: /*
79: Set function evaluation routine and vector
80: */
81: SNESSetFunction(snes,r,FormFunction,(void*)F);
84: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85: Create matrix data structure; set Jacobian evaluation routine
86: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
88: MatCreate(PETSC_COMM_WORLD,&J);
89: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);
90: MatSetFromOptions(J);
92: /*
93: Set Jacobian matrix data structure and default Jacobian evaluation
94: routine. User can override with:
95: -snes_fd : default finite differencing approximation of Jacobian
96: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
97: (unless user explicitly sets preconditioner)
98: -snes_mf_operator : form preconditioning matrix as set by the user,
99: but use matrix-free approx for Jacobian-vector
100: products within Newton-Krylov method
101: */
103: SNESSetJacobian(snes,J,J,FormJacobian,PETSC_NULL);
105: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106: Customize nonlinear solver; set runtime options
107: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109: /*
110: Set an optional user-defined monitoring routine
111: */
112: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,400,400,&monP.viewer);
113: SNESMonitorSet(snes,Monitor,&monP,0);
115: /*
116: Set names for some vectors to facilitate monitoring (optional)
117: */
118: PetscObjectSetName((PetscObject)x,"Approximate Solution");
119: PetscObjectSetName((PetscObject)U,"Exact Solution");
121: /*
122: Set SNES/KSP/KSP/PC runtime options, e.g.,
123: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
124: */
125: SNESSetFromOptions(snes);
127: /*
128: Print parameters used for convergence testing (optional) ... just
129: to demonstrate this routine; this information is also printed with
130: the option -snes_view
131: */
132: SNESGetTolerances(snes,&abstol,&rtol,&stol,&maxit,&maxf);
133: PetscPrintf(PETSC_COMM_WORLD,"atol=%G, rtol=%G, stol=%G, maxit=%D, maxf=%D\n",abstol,rtol,stol,maxit,maxf);
135: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136: Initialize application:
137: Store right-hand-side of PDE and exact solution
138: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
140: xp = 0.0;
141: for (i=0; i<n; i++) {
142: v = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
143: VecSetValues(F,1,&i,&v,INSERT_VALUES);
144: v= xp*xp*xp;
145: VecSetValues(U,1,&i,&v,INSERT_VALUES);
146: xp += h;
147: }
149: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: Evaluate initial guess; then solve nonlinear system
151: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
152: /*
153: Note: The user should initialize the vector, x, with the initial guess
154: for the nonlinear solver prior to calling SNESSolve(). In particular,
155: to employ an initial guess of zero, the user should explicitly set
156: this vector to zero by calling VecSet().
157: */
158: FormInitialGuess(x);
159: SNESSolve(snes,PETSC_NULL,x);
160: SNESGetIterationNumber(snes,&its);
161: PetscPrintf(PETSC_COMM_WORLD,"number of Newton iterations = %D\n\n",its);
163: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164: Check solution and clean up
165: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
167: /*
168: Check the error
169: */
170: VecAXPY(x,none,U);
171: VecNorm(x,NORM_2,&norm);
172: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %A, Iterations %D\n",norm,its);
175: /*
176: Free work space. All PETSc objects should be destroyed when they
177: are no longer needed.
178: */
179: VecDestroy(x); VecDestroy(r);
180: VecDestroy(U); VecDestroy(F);
181: MatDestroy(J); SNESDestroy(snes);
182: PetscViewerDestroy(monP.viewer);
183: PetscFinalize();
185: return 0;
186: }
187: /* ------------------------------------------------------------------- */
190: /*
191: FormInitialGuess - Computes initial guess.
193: Input/Output Parameter:
194: . x - the solution vector
195: */
196: PetscErrorCode FormInitialGuess(Vec x)
197: {
199: PetscScalar pfive = .50;
200: VecSet(x,pfive);
201: return 0;
202: }
203: /* ------------------------------------------------------------------- */
206: /*
207: FormFunction - Evaluates nonlinear function, F(x).
209: Input Parameters:
210: . snes - the SNES context
211: . x - input vector
212: . ctx - optional user-defined context, as set by SNESSetFunction()
214: Output Parameter:
215: . f - function vector
217: Note:
218: The user-defined context can contain any application-specific data
219: needed for the function evaluation (such as various parameters, work
220: vectors, and grid information). In this program the context is just
221: a vector containing the right-hand-side of the discretized PDE.
222: */
224: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
225: {
226: Vec g = (Vec)ctx;
227: PetscScalar *xx,*ff,*gg,d;
229: PetscInt i,n;
231: /*
232: Get pointers to vector data.
233: - For default PETSc vectors, VecGetArray() returns a pointer to
234: the data array. Otherwise, the routine is implementation dependent.
235: - You MUST call VecRestoreArray() when you no longer need access to
236: the array.
237: */
238: VecGetArray(x,&xx);
239: VecGetArray(f,&ff);
240: VecGetArray(g,&gg);
242: /*
243: Compute function
244: */
245: VecGetSize(x,&n);
246: d = (PetscReal)(n - 1); d = d*d;
247: ff[0] = xx[0];
248: for (i=1; i<n-1; i++) {
249: ff[i] = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - gg[i];
250: }
251: ff[n-1] = xx[n-1] - 1.0;
253: /*
254: Restore vectors
255: */
256: VecRestoreArray(x,&xx);
257: VecRestoreArray(f,&ff);
258: VecRestoreArray(g,&gg);
259: return 0;
260: }
261: /* ------------------------------------------------------------------- */
264: /*
265: FormJacobian - Evaluates Jacobian matrix.
267: Input Parameters:
268: . snes - the SNES context
269: . x - input vector
270: . dummy - optional user-defined context (not used here)
272: Output Parameters:
273: . jac - Jacobian matrix
274: . B - optionally different preconditioning matrix
275: . flag - flag indicating matrix structure
276: */
278: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure*flag,void *dummy)
279: {
280: PetscScalar *xx,A[3],d;
282: PetscInt i,n,j[3];
284: /*
285: Get pointer to vector data
286: */
287: VecGetArray(x,&xx);
289: /*
290: Compute Jacobian entries and insert into matrix.
291: - Note that in this case we set all elements for a particular
292: row at once.
293: */
294: VecGetSize(x,&n);
295: d = (PetscReal)(n - 1); d = d*d;
297: /*
298: Interior grid points
299: */
300: for (i=1; i<n-1; i++) {
301: j[0] = i - 1; j[1] = i; j[2] = i + 1;
302: A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
303: MatSetValues(*jac,1,&i,3,j,A,INSERT_VALUES);
304: }
306: /*
307: Boundary points
308: */
309: i = 0; A[0] = 1.0;
310: MatSetValues(*jac,1,&i,1,&i,A,INSERT_VALUES);
311: i = n-1; A[0] = 1.0;
312: MatSetValues(*jac,1,&i,1,&i,A,INSERT_VALUES);
314: /*
315: Restore vector
316: */
317: VecRestoreArray(x,&xx);
319: /*
320: Assemble matrix
321: */
322: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
323: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
325: return 0;
326: }
327: /* ------------------------------------------------------------------- */
330: /*
331: Monitor - User-defined monitoring routine that views the
332: current iterate with an x-window plot.
334: Input Parameters:
335: snes - the SNES context
336: its - iteration number
337: norm - 2-norm function value (may be estimated)
338: ctx - optional user-defined context for private data for the
339: monitor routine, as set by SNESMonitorSet()
341: Note:
342: See the manpage for PetscViewerDrawOpen() for useful runtime options,
343: such as -nox to deactivate all x-window output.
344: */
345: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
346: {
348: MonitorCtx *monP = (MonitorCtx*) ctx;
349: Vec x;
351: PetscPrintf(PETSC_COMM_WORLD,"iter = %D, SNES Function norm %G\n",its,fnorm);
352: SNESGetSolution(snes,&x);
353: VecView(x,monP->viewer);
354: return 0;
355: }