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.nn";
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: petscsles.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;
40: #undef __FUNCT__
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/SLES/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=%dn",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,&its);
158: PetscPrintf(PETSC_COMM_WORLD,"number of Newton iterations = %dnn",its);
160: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
161: Check solution and clean up
162: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
164: /*
165: Check the error
166: */
167: VecAXPY(&none,U,x);
168: ierr = VecNorm(x,NORM_2,&norm);
169: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %A, Iterations %dn",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();
182: return 0;
183: }
184: /* ------------------------------------------------------------------- */
185: #undef __FUNCT__
187: /*
188: FormInitialGuess - Computes initial guess.
190: Input/Output Parameter:
191: . x - the solution vector
192: */
193: int FormInitialGuess(Vec x)
194: {
195: int ierr;
196: PetscScalar pfive = .50;
197: VecSet(&pfive,x);
198: return 0;
199: }
200: /* ------------------------------------------------------------------- */
201: #undef __FUNCT__
203: /*
204: FormFunction - Evaluates nonlinear function, F(x).
206: Input Parameters:
207: . snes - the SNES context
208: . x - input vector
209: . ctx - optional user-defined context, as set by SNESSetFunction()
211: Output Parameter:
212: . f - function vector
214: Note:
215: The user-defined context can contain any application-specific data
216: needed for the function evaluation (such as various parameters, work
217: vectors, and grid information). In this program the context is just
218: a vector containing the right-hand-side of the discretized PDE.
219: */
221: int FormFunction(SNES snes,Vec x,Vec f,void *ctx)
222: {
223: Vec g = (Vec)ctx;
224: PetscScalar *xx,*ff,*gg,d;
225: int i,ierr,n;
227: /*
228: Get pointers to vector data.
229: - For default PETSc vectors, VecGetArray() returns a pointer to
230: the data array. Otherwise, the routine is implementation dependent.
231: - You MUST call VecRestoreArray() when you no longer need access to
232: the array.
233: */
234: VecGetArray(x,&xx);
235: VecGetArray(f,&ff);
236: VecGetArray(g,&gg);
238: /*
239: Compute function
240: */
241: VecGetSize(x,&n);
242: d = (PetscReal)(n - 1); d = d*d;
243: ff[0] = xx[0];
244: for (i=1; i<n-1; i++) {
245: ff[i] = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - gg[i];
246: }
247: ff[n-1] = xx[n-1] - 1.0;
249: /*
250: Restore vectors
251: */
252: VecRestoreArray(x,&xx);
253: VecRestoreArray(f,&ff);
254: VecRestoreArray(g,&gg);
255: return 0;
256: }
257: /* ------------------------------------------------------------------- */
258: #undef __FUNCT__
260: /*
261: FormJacobian - Evaluates Jacobian matrix.
263: Input Parameters:
264: . snes - the SNES context
265: . x - input vector
266: . dummy - optional user-defined context (not used here)
268: Output Parameters:
269: . jac - Jacobian matrix
270: . B - optionally different preconditioning matrix
271: . flag - flag indicating matrix structure
272: */
274: int FormJacobian(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure*flag,void *dummy)
275: {
276: PetscScalar *xx,A[3],d;
277: int i,n,j[3],ierr;
279: /*
280: Get pointer to vector data
281: */
282: VecGetArray(x,&xx);
284: /*
285: Compute Jacobian entries and insert into matrix.
286: - Note that in this case we set all elements for a particular
287: row at once.
288: */
289: VecGetSize(x,&n);
290: d = (PetscReal)(n - 1); d = d*d;
292: /*
293: Interior grid points
294: */
295: for (i=1; i<n-1; i++) {
296: j[0] = i - 1; j[1] = i; j[2] = i + 1;
297: A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
298: MatSetValues(*jac,1,&i,3,j,A,INSERT_VALUES);
299: }
301: /*
302: Boundary points
303: */
304: i = 0; A[0] = 1.0;
305: MatSetValues(*jac,1,&i,1,&i,A,INSERT_VALUES);
306: i = n-1; A[0] = 1.0;
307: MatSetValues(*jac,1,&i,1,&i,A,INSERT_VALUES);
309: /*
310: Restore vector
311: */
312: VecRestoreArray(x,&xx);
314: /*
315: Assemble matrix
316: */
317: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
318: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
320: return 0;
321: }
322: /* ------------------------------------------------------------------- */
323: #undef __FUNCT__
325: /*
326: Monitor - User-defined monitoring routine that views the
327: current iterate with an x-window plot.
329: Input Parameters:
330: snes - the SNES context
331: its - iteration number
332: norm - 2-norm function value (may be estimated)
333: ctx - optional user-defined context for private data for the
334: monitor routine, as set by SNESSetMonitor()
336: Note:
337: See the manpage for PetscViewerDrawOpen() for useful runtime options,
338: such as -nox to deactivate all x-window output.
339: */
340: int Monitor(SNES snes,int its,PetscReal fnorm,void *ctx)
341: {
342: int ierr;
343: MonitorCtx *monP = (MonitorCtx*) ctx;
344: Vec x;
346: PetscPrintf(PETSC_COMM_WORLD,"iter = %d, SNES Function norm %gn",its,fnorm);
347: SNESGetSolution(snes,&x);
348: VecView(x,monP->viewer);
349: return 0;
350: }