Actual source code: ex2.c
1: /*$Id: ex2.c,v 1.83 2001/08/07 21:31:17 bsmith Exp $*/
3: static char help[] = "Newton method to solve u'' + u^{2} = f, sequentially.\n\
4: This example employs a user-defined monitoring routine.\n\n";
6: /*T
7: Concepts: SNES^basic uniprocessor example
8: Concepts: SNES^setting a user-defined monitoring routine
9: Processors: 1
10: T*/
12: /*
13: Include "petscdraw.h" so that we can use PETSc drawing routines.
14: Include "petscsnes.h" so that we can use SNES solvers. Note that this
15: file automatically includes:
16: petsc.h - base PETSc routines petscvec.h - vectors
17: petscsys.h - system routines petscmat.h - matrices
18: petscis.h - index sets petscksp.h - Krylov subspace methods
19: petscviewer.h - viewers petscpc.h - preconditioners
20: petscksp.h - linear solvers
21: */
23: #include petscsnes.h
25: /*
26: User-defined routines
27: */
28: extern int FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
29: extern int FormFunction(SNES,Vec,Vec,void*);
30: extern int FormInitialGuess(Vec);
31: extern int Monitor(SNES,int,PetscReal,void *);
33: /*
34: User-defined context for monitoring
35: */
36: typedef struct {
37: PetscViewer viewer;
38: } MonitorCtx;
42: int main(int argc,char **argv)
43: {
44: SNES snes; /* SNES context */
45: Vec x,r,F,U; /* vectors */
46: Mat J; /* Jacobian matrix */
47: MonitorCtx monP; /* monitoring context */
48: int ierr,its,n = 5,i,maxit,maxf,size;
49: PetscScalar h,xp,v,none = -1.0;
50: PetscReal atol,rtol,stol,norm;
52: PetscInitialize(&argc,&argv,(char *)0,help);
53: MPI_Comm_size(PETSC_COMM_WORLD,&size);
54: if (size != 1) SETERRQ(1,"This is a uniprocessor example only!");
55: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
56: h = 1.0/(n-1);
58: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
59: Create nonlinear solver context
60: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
62: SNESCreate(PETSC_COMM_WORLD,&snes);
64: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
65: Create vector data structures; set function evaluation routine
66: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
67: /*
68: Note that we form 1 vector from scratch and then duplicate as needed.
69: */
70: VecCreate(PETSC_COMM_WORLD,&x);
71: VecSetSizes(x,PETSC_DECIDE,n);
72: VecSetFromOptions(x);
73: VecDuplicate(x,&r);
74: VecDuplicate(x,&F);
75: VecDuplicate(x,&U);
77: /*
78: Set function evaluation routine and vector
79: */
80: SNESSetFunction(snes,r,FormFunction,(void*)F);
83: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84: Create matrix data structure; set Jacobian evaluation routine
85: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
87: MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,&J);
88: MatSetFromOptions(J);
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,PETSC_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: SNESSetMonitor(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,&atol,&rtol,&stol,&maxit,&maxf);
131: PetscPrintf(PETSC_COMM_WORLD,"atol=%g, rtol=%g, stol=%g, maxit=%d, maxf=%d\n",atol,rtol,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,x);
158: SNESGetIterationNumber(snes,&its);
159: PetscPrintf(PETSC_COMM_WORLD,"number of Newton iterations = %d\n\n",its);
161: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162: Check solution and clean up
163: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
165: /*
166: Check the error
167: */
168: VecAXPY(&none,U,x);
169: VecNorm(x,NORM_2,&norm);
170: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %A, Iterations %d\n",norm,its);
173: /*
174: Free work space. All PETSc objects should be destroyed when they
175: are no longer needed.
176: */
177: VecDestroy(x); VecDestroy(r);
178: VecDestroy(U); VecDestroy(F);
179: MatDestroy(J); SNESDestroy(snes);
180: PetscViewerDestroy(monP.viewer);
181: PetscFinalize();
183: return 0;
184: }
185: /* ------------------------------------------------------------------- */
188: /*
189: FormInitialGuess - Computes initial guess.
191: Input/Output Parameter:
192: . x - the solution vector
193: */
194: int FormInitialGuess(Vec x)
195: {
196: int ierr;
197: PetscScalar pfive = .50;
198: VecSet(&pfive,x);
199: return 0;
200: }
201: /* ------------------------------------------------------------------- */
204: /*
205: FormFunction - Evaluates nonlinear function, F(x).
207: Input Parameters:
208: . snes - the SNES context
209: . x - input vector
210: . ctx - optional user-defined context, as set by SNESSetFunction()
212: Output Parameter:
213: . f - function vector
215: Note:
216: The user-defined context can contain any application-specific data
217: needed for the function evaluation (such as various parameters, work
218: vectors, and grid information). In this program the context is just
219: a vector containing the right-hand-side of the discretized PDE.
220: */
222: int FormFunction(SNES snes,Vec x,Vec f,void *ctx)
223: {
224: Vec g = (Vec)ctx;
225: PetscScalar *xx,*ff,*gg,d;
226: int i,ierr,n;
228: /*
229: Get pointers to vector data.
230: - For default PETSc vectors, VecGetArray() returns a pointer to
231: the data array. Otherwise, the routine is implementation dependent.
232: - You MUST call VecRestoreArray() when you no longer need access to
233: the array.
234: */
235: VecGetArray(x,&xx);
236: VecGetArray(f,&ff);
237: VecGetArray(g,&gg);
239: /*
240: Compute function
241: */
242: VecGetSize(x,&n);
243: d = (PetscReal)(n - 1); d = d*d;
244: ff[0] = xx[0];
245: for (i=1; i<n-1; i++) {
246: ff[i] = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - gg[i];
247: }
248: ff[n-1] = xx[n-1] - 1.0;
250: /*
251: Restore vectors
252: */
253: VecRestoreArray(x,&xx);
254: VecRestoreArray(f,&ff);
255: VecRestoreArray(g,&gg);
256: return 0;
257: }
258: /* ------------------------------------------------------------------- */
261: /*
262: FormJacobian - Evaluates Jacobian matrix.
264: Input Parameters:
265: . snes - the SNES context
266: . x - input vector
267: . dummy - optional user-defined context (not used here)
269: Output Parameters:
270: . jac - Jacobian matrix
271: . B - optionally different preconditioning matrix
272: . flag - flag indicating matrix structure
273: */
275: int FormJacobian(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure*flag,void *dummy)
276: {
277: PetscScalar *xx,A[3],d;
278: int i,n,j[3],ierr;
280: /*
281: Get pointer to vector data
282: */
283: VecGetArray(x,&xx);
285: /*
286: Compute Jacobian entries and insert into matrix.
287: - Note that in this case we set all elements for a particular
288: row at once.
289: */
290: VecGetSize(x,&n);
291: d = (PetscReal)(n - 1); d = d*d;
293: /*
294: Interior grid points
295: */
296: for (i=1; i<n-1; i++) {
297: j[0] = i - 1; j[1] = i; j[2] = i + 1;
298: A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
299: MatSetValues(*jac,1,&i,3,j,A,INSERT_VALUES);
300: }
302: /*
303: Boundary points
304: */
305: i = 0; A[0] = 1.0;
306: MatSetValues(*jac,1,&i,1,&i,A,INSERT_VALUES);
307: i = n-1; A[0] = 1.0;
308: MatSetValues(*jac,1,&i,1,&i,A,INSERT_VALUES);
310: /*
311: Restore vector
312: */
313: VecRestoreArray(x,&xx);
315: /*
316: Assemble matrix
317: */
318: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
319: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
321: return 0;
322: }
323: /* ------------------------------------------------------------------- */
326: /*
327: Monitor - User-defined monitoring routine that views the
328: current iterate with an x-window plot.
330: Input Parameters:
331: snes - the SNES context
332: its - iteration number
333: norm - 2-norm function value (may be estimated)
334: ctx - optional user-defined context for private data for the
335: monitor routine, as set by SNESSetMonitor()
337: Note:
338: See the manpage for PetscViewerDrawOpen() for useful runtime options,
339: such as -nox to deactivate all x-window output.
340: */
341: int Monitor(SNES snes,int its,PetscReal fnorm,void *ctx)
342: {
343: int ierr;
344: MonitorCtx *monP = (MonitorCtx*) ctx;
345: Vec x;
347: PetscPrintf(PETSC_COMM_WORLD,"iter = %d, SNES Function norm %g\n",its,fnorm);
348: SNESGetSolution(snes,&x);
349: VecView(x,monP->viewer);
350: return 0;
351: }