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