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 */
46: PetscInt its,n = 5,i,maxit,maxf;
47: PetscMPIInt size;
48: PetscScalar h,xp,v,none = -1.0;
49: PetscReal abstol,rtol,stol,norm;
51: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
52: MPI_Comm_size(PETSC_COMM_WORLD,&size);
53: if (size != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This is a uniprocessor example only!");
54: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
55: h = 1.0/(n-1);
57: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
58: Create nonlinear solver context
59: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
61: SNESCreate(PETSC_COMM_WORLD,&snes);
63: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
64: Create vector data structures; set function evaluation routine
65: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
66: /*
67: Note that we form 1 vector from scratch and then duplicate as needed.
68: */
69: VecCreate(PETSC_COMM_WORLD,&x);
70: VecSetSizes(x,PETSC_DECIDE,n);
71: VecSetFromOptions(x);
72: VecDuplicate(x,&r);
73: VecDuplicate(x,&F);
74: VecDuplicate(x,&U);
76: /*
77: Set function evaluation routine and vector
78: */
79: SNESSetFunction(snes,r,FormFunction,(void*)F);
81: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
82: Create matrix data structure; set Jacobian evaluation routine
83: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
85: MatCreate(PETSC_COMM_WORLD,&J);
86: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);
87: MatSetFromOptions(J);
88: MatSeqAIJSetPreallocation(J,3,NULL);
90: /*
91: Set Jacobian matrix data structure and default Jacobian evaluation
92: routine. User can override with:
93: -snes_fd : default finite differencing approximation of Jacobian
94: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
95: (unless user explicitly sets preconditioner)
96: -snes_mf_operator : form preconditioning matrix as set by the user,
97: but use matrix-free approx for Jacobian-vector
98: products within Newton-Krylov method
99: */
101: SNESSetJacobian(snes,J,J,FormJacobian,NULL);
103: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104: Customize nonlinear solver; set runtime options
105: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
107: /*
108: Set an optional user-defined monitoring routine
109: */
110: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,400,400,&monP.viewer);
111: SNESMonitorSet(snes,Monitor,&monP,0);
113: /*
114: Set names for some vectors to facilitate monitoring (optional)
115: */
116: PetscObjectSetName((PetscObject)x,"Approximate Solution");
117: PetscObjectSetName((PetscObject)U,"Exact Solution");
119: /*
120: Set SNES/KSP/KSP/PC runtime options, e.g.,
121: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
122: */
123: SNESSetFromOptions(snes);
125: /*
126: Print parameters used for convergence testing (optional) ... just
127: to demonstrate this routine; this information is also printed with
128: the option -snes_view
129: */
130: SNESGetTolerances(snes,&abstol,&rtol,&stol,&maxit,&maxf);
131: PetscPrintf(PETSC_COMM_WORLD,"atol=%g, rtol=%g, stol=%g, maxit=%D, maxf=%D\n",(double)abstol,(double)rtol,(double)stol,maxit,maxf);
133: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134: Initialize application:
135: Store right-hand-side of PDE and exact solution
136: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
138: xp = 0.0;
139: for (i=0; i<n; i++) {
140: v = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
141: VecSetValues(F,1,&i,&v,INSERT_VALUES);
142: v = xp*xp*xp;
143: VecSetValues(U,1,&i,&v,INSERT_VALUES);
144: xp += h;
145: }
147: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148: Evaluate initial guess; then solve nonlinear system
149: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
150: /*
151: Note: The user should initialize the vector, x, with the initial guess
152: for the nonlinear solver prior to calling SNESSolve(). In particular,
153: to employ an initial guess of zero, the user should explicitly set
154: this vector to zero by calling VecSet().
155: */
156: FormInitialGuess(x);
157: SNESSolve(snes,NULL,x);
158: SNESGetIterationNumber(snes,&its);
159: PetscPrintf(PETSC_COMM_WORLD,"number of SNES iterations = %D\n\n",its);
161: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162: Check solution and clean up
163: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
165: /*
166: Check the error
167: */
168: VecAXPY(x,none,U);
169: VecNorm(x,NORM_2,&norm);
170: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);
172: /*
173: Free work space. All PETSc objects should be destroyed when they
174: are no longer needed.
175: */
176: VecDestroy(&x); VecDestroy(&r);
177: VecDestroy(&U); VecDestroy(&F);
178: MatDestroy(&J); SNESDestroy(&snes);
179: PetscViewerDestroy(&monP.viewer);
180: PetscFinalize();
181: return ierr;
182: }
183: /* ------------------------------------------------------------------- */
184: /*
185: FormInitialGuess - Computes initial guess.
187: Input/Output Parameter:
188: . x - the solution vector
189: */
190: PetscErrorCode FormInitialGuess(Vec x)
191: {
193: PetscScalar pfive = .50;
194: VecSet(x,pfive);
195: return 0;
196: }
197: /* ------------------------------------------------------------------- */
198: /*
199: FormFunction - Evaluates nonlinear function, F(x).
201: Input Parameters:
202: . snes - the SNES context
203: . x - input vector
204: . ctx - optional user-defined context, as set by SNESSetFunction()
206: Output Parameter:
207: . f - function vector
209: Note:
210: The user-defined context can contain any application-specific data
211: needed for the function evaluation (such as various parameters, work
212: vectors, and grid information). In this program the context is just
213: a vector containing the right-hand-side of the discretized PDE.
214: */
216: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
217: {
218: Vec g = (Vec)ctx;
219: const PetscScalar *xx,*gg;
220: PetscScalar *ff,d;
221: PetscErrorCode ierr;
222: PetscInt i,n;
224: /*
225: Get pointers to vector data.
226: - For default PETSc vectors, VecGetArray() returns a pointer to
227: the data array. Otherwise, the routine is implementation dependent.
228: - You MUST call VecRestoreArray() when you no longer need access to
229: the array.
230: */
231: VecGetArrayRead(x,&xx);
232: VecGetArray(f,&ff);
233: VecGetArrayRead(g,&gg);
235: /*
236: Compute function
237: */
238: VecGetSize(x,&n);
239: d = (PetscReal)(n - 1); d = d*d;
240: ff[0] = xx[0];
241: 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];
242: ff[n-1] = xx[n-1] - 1.0;
244: /*
245: Restore vectors
246: */
247: VecRestoreArrayRead(x,&xx);
248: VecRestoreArray(f,&ff);
249: VecRestoreArrayRead(g,&gg);
250: return 0;
251: }
252: /* ------------------------------------------------------------------- */
253: /*
254: FormJacobian - Evaluates Jacobian matrix.
256: Input Parameters:
257: . snes - the SNES context
258: . x - input vector
259: . dummy - optional user-defined context (not used here)
261: Output Parameters:
262: . jac - Jacobian matrix
263: . B - optionally different preconditioning matrix
265: */
267: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
268: {
269: const PetscScalar *xx;
270: PetscScalar A[3],d;
271: PetscErrorCode ierr;
272: PetscInt i,n,j[3];
274: /*
275: Get pointer to vector data
276: */
277: VecGetArrayRead(x,&xx);
279: /*
280: Compute Jacobian entries and insert into matrix.
281: - Note that in this case we set all elements for a particular
282: row at once.
283: */
284: VecGetSize(x,&n);
285: d = (PetscReal)(n - 1); d = d*d;
287: /*
288: Interior grid points
289: */
290: for (i=1; i<n-1; i++) {
291: j[0] = i - 1; j[1] = i; j[2] = i + 1;
292: A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
293: MatSetValues(B,1,&i,3,j,A,INSERT_VALUES);
294: }
296: /*
297: Boundary points
298: */
299: i = 0; A[0] = 1.0;
301: MatSetValues(B,1,&i,1,&i,A,INSERT_VALUES);
303: i = n-1; A[0] = 1.0;
305: MatSetValues(B,1,&i,1,&i,A,INSERT_VALUES);
307: /*
308: Restore vector
309: */
310: VecRestoreArrayRead(x,&xx);
312: /*
313: Assemble matrix
314: */
315: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
316: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
317: if (jac != B) {
318: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
319: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
320: }
321: return 0;
322: }
323: /* ------------------------------------------------------------------- */
324: /*
325: Monitor - User-defined monitoring routine that views the
326: current iterate with an x-window plot.
328: Input Parameters:
329: snes - the SNES context
330: its - iteration number
331: norm - 2-norm function value (may be estimated)
332: ctx - optional user-defined context for private data for the
333: monitor routine, as set by SNESMonitorSet()
335: Note:
336: See the manpage for PetscViewerDrawOpen() for useful runtime options,
337: such as -nox to deactivate all x-window output.
338: */
339: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
340: {
342: MonitorCtx *monP = (MonitorCtx*) ctx;
343: Vec x;
345: PetscPrintf(PETSC_COMM_WORLD,"iter = %D, SNES Function norm %g\n",its,(double)fnorm);
346: SNESGetSolution(snes,&x);
347: VecView(x,monP->viewer);
348: return 0;
349: }
351: /*TEST
353: test:
354: args: -nox -snes_monitor_cancel -snes_monitor_short -snes_view -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always
356: test:
357: suffix: 2
358: args: -nox -snes_monitor_cancel -snes_monitor_short -snes_type newtontr -snes_view
359: requires: !single
361: test:
362: suffix: 3
363: 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
365: TEST*/