Actual source code: ex1.c
1: /*$Id: ex1.c,v 1.87 2001/08/27 17:31:28 bsmith Exp $*/
3: /* Program usage: ex4 [-help] [all PETSc options] */
5: static char help[] = "Solves a nonlinear system on 1 processor with SNES. We\n\
6: solve the Bratu (SFI - solid fuel ignition) problem in a 2D rectangular domain.\n\
7: This example also illustrates the use of matrix coloring. Runtime options include:\n\
8: -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
9: problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81)\n\
10: -mx <xg>, where <xg> = number of grid points in the x-direction\n\
11: -my <yg>, where <yg> = number of grid points in the y-direction\n\n";
13: /*T
14: Concepts: SNES^sequential Bratu example
15: Processors: 1
16: T*/
18: /* ------------------------------------------------------------------------
20: Solid Fuel Ignition (SFI) problem. This problem is modeled by
21: the partial differential equation
22:
23: -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
24:
25: with boundary conditions
26:
27: u = 0 for x = 0, x = 1, y = 0, y = 1.
28:
29: A finite difference approximation with the usual 5-point stencil
30: is used to discretize the boundary value problem to obtain a nonlinear
31: system of equations.
33: The parallel version of this code is snes/examples/tutorials/ex5.c
35: ------------------------------------------------------------------------- */
37: /*
38: Include "petscsnes.h" so that we can use SNES solvers. Note that
39: this file automatically includes:
40: petsc.h - base PETSc routines petscvec.h - vectors
41: petscsys.h - system routines petscmat.h - matrices
42: petscis.h - index sets petscksp.h - Krylov subspace methods
43: petscviewer.h - viewers petscpc.h - preconditioners
44: petscksp.h - linear solvers
45: */
47: #include petscsnes.h
49: /*
50: User-defined application context - contains data needed by the
51: application-provided call-back routines, FormJacobian() and
52: FormFunction().
53: */
54: typedef struct {
55: PetscReal param; /* test problem parameter */
56: int mx; /* Discretization in x-direction */
57: int my; /* Discretization in y-direction */
58: } AppCtx;
60: /*
61: User-defined routines
62: */
63: extern int FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
64: extern int FormFunction(SNES,Vec,Vec,void*);
65: extern int FormInitialGuess(AppCtx*,Vec);
69: int main(int argc,char **argv)
70: {
71: SNES snes; /* nonlinear solver context */
72: Vec x,r; /* solution, residual vectors */
73: Mat J; /* Jacobian matrix */
74: AppCtx user; /* user-defined application context */
75: int i,ierr,its,N,size,hist_its[50];
76: PetscReal bratu_lambda_max = 6.81,bratu_lambda_min = 0.,history[50];
77: MatFDColoring fdcoloring;
78: PetscTruth matrix_free,flg,fd_coloring;
80: PetscInitialize(&argc,&argv,(char *)0,help);
81: MPI_Comm_size(PETSC_COMM_WORLD,&size);
82: if (size != 1) SETERRQ(1,"This is a uniprocessor example only!");
84: /*
85: Initialize problem parameters
86: */
87: user.mx = 4; user.my = 4; user.param = 6.0;
88: PetscOptionsGetInt(PETSC_NULL,"-mx",&user.mx,PETSC_NULL);
89: PetscOptionsGetInt(PETSC_NULL,"-my",&user.my,PETSC_NULL);
90: PetscOptionsGetReal(PETSC_NULL,"-par",&user.param,PETSC_NULL);
91: if (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min) {
92: SETERRQ(1,"Lambda is out of range");
93: }
94: N = user.mx*user.my;
95:
96: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97: Create nonlinear solver context
98: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
100: SNESCreate(PETSC_COMM_WORLD,&snes);
102: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103: Create vector data structures; set function evaluation routine
104: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
106: VecCreate(PETSC_COMM_WORLD,&x);
107: VecSetSizes(x,PETSC_DECIDE,N);
108: VecSetFromOptions(x);
109: VecDuplicate(x,&r);
111: /*
112: Set function evaluation routine and vector. Whenever the nonlinear
113: solver needs to evaluate the nonlinear function, it will call this
114: routine.
115: - Note that the final routine argument is the user-defined
116: context that provides application-specific data for the
117: function evaluation routine.
118: */
119: SNESSetFunction(snes,r,FormFunction,(void *)&user);
121: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122: Create matrix data structure; set Jacobian evaluation routine
123: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
125: /*
126: Create matrix. Here we only approximately preallocate storage space
127: for the Jacobian. See the users manual for a discussion of better
128: techniques for preallocating matrix memory.
129: */
130: PetscOptionsHasName(PETSC_NULL,"-snes_mf",&matrix_free);
131: if (!matrix_free) {
132: MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,5,PETSC_NULL,&J);
133: }
135: /*
136: This option will cause the Jacobian to be computed via finite differences
137: efficiently using a coloring of the columns of the matrix.
138: */
139: PetscOptionsHasName(PETSC_NULL,"-snes_fd_coloring",&fd_coloring);
141: if (matrix_free && fd_coloring) SETERRQ(1,"Use only one of -snes_mf, -snes_fd_coloring options!\n\
142: You can do -snes_mf_operator -snes_fd_coloring");
144: if (fd_coloring) {
145: ISColoring iscoloring;
146: MatStructure str;
148: /*
149: This initializes the nonzero structure of the Jacobian. This is artificial
150: because clearly if we had a routine to compute the Jacobian we won't need
151: to use finite differences.
152: */
153: FormJacobian(snes,x,&J,&J,&str,&user);
155: /*
156: Color the matrix, i.e. determine groups of columns that share no common
157: rows. These columns in the Jacobian can all be computed simulataneously.
158: */
159: MatGetColoring(J,MATCOLORING_NATURAL,&iscoloring);
160: /*
161: Create the data structure that SNESDefaultComputeJacobianColor() uses
162: to compute the actual Jacobians via finite differences.
163: */
164: MatFDColoringCreate(J,iscoloring,&fdcoloring);
165: MatFDColoringSetFunction(fdcoloring,(int (*)(void))FormFunction,&user);
166: MatFDColoringSetFromOptions(fdcoloring);
167: /*
168: Tell SNES to use the routine SNESDefaultComputeJacobianColor()
169: to compute Jacobians.
170: */
171: SNESSetJacobian(snes,J,J,SNESDefaultComputeJacobianColor,fdcoloring);
172: ISColoringDestroy(iscoloring);
173: }
174: /*
175: Set Jacobian matrix data structure and default Jacobian evaluation
176: routine. Whenever the nonlinear solver needs to compute the
177: Jacobian matrix, it will call this routine.
178: - Note that the final routine argument is the user-defined
179: context that provides application-specific data for the
180: Jacobian evaluation routine.
181: - The user can override with:
182: -snes_fd : default finite differencing approximation of Jacobian
183: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
184: (unless user explicitly sets preconditioner)
185: -snes_mf_operator : form preconditioning matrix as set by the user,
186: but use matrix-free approx for Jacobian-vector
187: products within Newton-Krylov method
188: */
189: else if (!matrix_free) {
190: SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);
191: }
193: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
194: Customize nonlinear solver; set runtime options
195: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
197: /*
198: Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
199: */
200: SNESSetFromOptions(snes);
202: /*
203: Set array that saves the function norms. This array is intended
204: when the user wants to save the convergence history for later use
205: rather than just to view the function norms via -snes_monitor.
206: */
207: SNESSetConvergenceHistory(snes,history,hist_its,50,PETSC_TRUE);
209: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
210: Evaluate initial guess; then solve nonlinear system
211: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
212: /*
213: Note: The user should initialize the vector, x, with the initial guess
214: for the nonlinear solver prior to calling SNESSolve(). In particular,
215: to employ an initial guess of zero, the user should explicitly set
216: this vector to zero by calling VecSet().
217: */
218: FormInitialGuess(&user,x);
219: SNESSolve(snes,x);
220: SNESGetIterationNumber(snes,&its);
221: PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %d\n",its);
224: /*
225: Print the convergence history. This is intended just to demonstrate
226: use of the data attained via SNESSetConvergenceHistory().
227: */
228: PetscOptionsHasName(PETSC_NULL,"-print_history",&flg);
229: if (flg) {
230: for (i=0; i<its+1; i++) {
231: PetscPrintf(PETSC_COMM_WORLD,"iteration %d: Linear iterations %d Function norm = %g\n",i,hist_its[i],history[i]);
232: }
233: }
235: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
236: Free work space. All PETSc objects should be destroyed when they
237: are no longer needed.
238: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
240: if (!matrix_free) {
241: MatDestroy(J);
242: }
243: if (fd_coloring) {
244: MatFDColoringDestroy(fdcoloring);
245: }
246: VecDestroy(x);
247: VecDestroy(r);
248: SNESDestroy(snes);
249: PetscFinalize();
251: return 0;
252: }
253: /* ------------------------------------------------------------------- */
256: /*
257: FormInitialGuess - Forms initial approximation.
259: Input Parameters:
260: user - user-defined application context
261: X - vector
263: Output Parameter:
264: X - vector
265: */
266: int FormInitialGuess(AppCtx *user,Vec X)
267: {
268: int i,j,row,mx,my,ierr;
269: PetscReal lambda,temp1,temp,hx,hy;
270: PetscScalar *x;
272: mx = user->mx;
273: my = user->my;
274: lambda = user->param;
276: hx = 1.0 / (PetscReal)(mx-1);
277: hy = 1.0 / (PetscReal)(my-1);
279: /*
280: Get a pointer to vector data.
281: - For default PETSc vectors, VecGetArray() returns a pointer to
282: the data array. Otherwise, the routine is implementation dependent.
283: - You MUST call VecRestoreArray() when you no longer need access to
284: the array.
285: */
286: VecGetArray(X,&x);
287: temp1 = lambda/(lambda + 1.0);
288: for (j=0; j<my; j++) {
289: temp = (PetscReal)(PetscMin(j,my-j-1))*hy;
290: for (i=0; i<mx; i++) {
291: row = i + j*mx;
292: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
293: x[row] = 0.0;
294: continue;
295: }
296: x[row] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));
297: }
298: }
300: /*
301: Restore vector
302: */
303: VecRestoreArray(X,&x);
304: return 0;
305: }
306: /* ------------------------------------------------------------------- */
309: /*
310: FormFunction - Evaluates nonlinear function, F(x).
312: Input Parameters:
313: . snes - the SNES context
314: . X - input vector
315: . ptr - optional user-defined context, as set by SNESSetFunction()
317: Output Parameter:
318: . F - function vector
319: */
320: int FormFunction(SNES snes,Vec X,Vec F,void *ptr)
321: {
322: AppCtx *user = (AppCtx*)ptr;
323: int ierr,i,j,row,mx,my;
324: PetscReal two = 2.0,one = 1.0,lambda,hx,hy,hxdhy,hydhx;
325: PetscScalar ut,ub,ul,ur,u,uxx,uyy,sc,*x,*f;
327: mx = user->mx;
328: my = user->my;
329: lambda = user->param;
330: hx = one / (PetscReal)(mx-1);
331: hy = one / (PetscReal)(my-1);
332: sc = hx*hy;
333: hxdhy = hx/hy;
334: hydhx = hy/hx;
336: /*
337: Get pointers to vector data
338: */
339: VecGetArray(X,&x);
340: VecGetArray(F,&f);
342: /*
343: Compute function
344: */
345: for (j=0; j<my; j++) {
346: for (i=0; i<mx; i++) {
347: row = i + j*mx;
348: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
349: f[row] = x[row];
350: continue;
351: }
352: u = x[row];
353: ub = x[row - mx];
354: ul = x[row - 1];
355: ut = x[row + mx];
356: ur = x[row + 1];
357: uxx = (-ur + two*u - ul)*hydhx;
358: uyy = (-ut + two*u - ub)*hxdhy;
359: f[row] = uxx + uyy - sc*lambda*PetscExpScalar(u);
360: }
361: }
363: /*
364: Restore vectors
365: */
366: VecRestoreArray(X,&x);
367: VecRestoreArray(F,&f);
368: return 0;
369: }
370: /* ------------------------------------------------------------------- */
373: /*
374: FormJacobian - Evaluates Jacobian matrix.
376: Input Parameters:
377: . snes - the SNES context
378: . x - input vector
379: . ptr - optional user-defined context, as set by SNESSetJacobian()
381: Output Parameters:
382: . A - Jacobian matrix
383: . B - optionally different preconditioning matrix
384: . flag - flag indicating matrix structure
385: */
386: int FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
387: {
388: AppCtx *user = (AppCtx*)ptr; /* user-defined applicatin context */
389: Mat jac = *J; /* Jacobian matrix */
390: int i,j,row,mx,my,col[5],ierr;
391: PetscScalar two = 2.0,one = 1.0,lambda,v[5],sc,*x;
392: PetscReal hx,hy,hxdhy,hydhx;
394: mx = user->mx;
395: my = user->my;
396: lambda = user->param;
397: hx = 1.0 / (PetscReal)(mx-1);
398: hy = 1.0 / (PetscReal)(my-1);
399: sc = hx*hy;
400: hxdhy = hx/hy;
401: hydhx = hy/hx;
403: /*
404: Get pointer to vector data
405: */
406: VecGetArray(X,&x);
408: /*
409: Compute entries of the Jacobian
410: */
411: for (j=0; j<my; j++) {
412: for (i=0; i<mx; i++) {
413: row = i + j*mx;
414: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
415: MatSetValues(jac,1,&row,1,&row,&one,INSERT_VALUES);
416: continue;
417: }
418: v[0] = -hxdhy; col[0] = row - mx;
419: v[1] = -hydhx; col[1] = row - 1;
420: v[2] = two*(hydhx + hxdhy) - sc*lambda*PetscExpScalar(x[row]); col[2] = row;
421: v[3] = -hydhx; col[3] = row + 1;
422: v[4] = -hxdhy; col[4] = row + mx;
423: MatSetValues(jac,1,&row,5,col,v,INSERT_VALUES);
424: }
425: }
427: /*
428: Restore vector
429: */
430: VecRestoreArray(X,&x);
432: /*
433: Assemble matrix
434: */
435: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
436: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
438: /*
439: Set flag to indicate that the Jacobian matrix retains an identical
440: nonzero structure throughout all nonlinear iterations (although the
441: values of the entries change). Thus, we can save some work in setting
442: up the preconditioner (e.g., no need to redo symbolic factorization for
443: ILU/ICC preconditioners).
444: - If the nonzero structure of the matrix is different during
445: successive linear solves, then the flag DIFFERENT_NONZERO_PATTERN
446: must be used instead. If you are unsure whether the matrix
447: structure has changed or not, use the flag DIFFERENT_NONZERO_PATTERN.
448: - Caution: If you specify SAME_NONZERO_PATTERN, PETSc
449: believes your assertion and does not check the structure
450: of the matrix. If you erroneously claim that the structure
451: is the same when it actually is not, the new preconditioner
452: will not function correctly. Thus, use this optimization
453: feature with caution!
454: */
455: *flag = SAME_NONZERO_PATTERN;
456: return 0;
457: }