Actual source code: ex1.c
2: static char help[] = "Solves the nonlinear system, the Bratu (SFI - solid fuel ignition) problem in a 2D rectangular domain.\n\
3: This example also illustrates the use of matrix coloring. Runtime options include:\n\
4: -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
5: problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81)\n\
6: -mx <xg>, where <xg> = number of grid points in the x-direction\n\
7: -my <yg>, where <yg> = number of grid points in the y-direction\n\n";
9: /*T
10: Concepts: SNES^sequential Bratu example
11: Processors: 1
12: T*/
14: /* ------------------------------------------------------------------------
16: Solid Fuel Ignition (SFI) problem. This problem is modeled by
17: the partial differential equation
19: -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
21: with boundary conditions
23: u = 0 for x = 0, x = 1, y = 0, y = 1.
25: A finite difference approximation with the usual 5-point stencil
26: is used to discretize the boundary value problem to obtain a nonlinear
27: system of equations.
29: The parallel version of this code is snes/tutorials/ex5.c
31: ------------------------------------------------------------------------- */
33: /*
34: Include "petscsnes.h" so that we can use SNES solvers. Note that
35: this file automatically includes:
36: petscsys.h - base PETSc routines petscvec.h - vectors
37: petscmat.h - matrices
38: petscis.h - index sets petscksp.h - Krylov subspace methods
39: petscviewer.h - viewers petscpc.h - preconditioners
40: petscksp.h - linear solvers
41: */
43: #include <petscsnes.h>
45: /*
46: User-defined application context - contains data needed by the
47: application-provided call-back routines, FormJacobian() and
48: FormFunction().
49: */
50: typedef struct {
51: PetscReal param; /* test problem parameter */
52: PetscInt mx; /* Discretization in x-direction */
53: PetscInt my; /* Discretization in y-direction */
54: } AppCtx;
56: /*
57: User-defined routines
58: */
59: extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
60: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
61: extern PetscErrorCode FormInitialGuess(AppCtx*,Vec);
62: extern PetscErrorCode ConvergenceTest(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*);
63: extern PetscErrorCode ConvergenceDestroy(void*);
64: extern PetscErrorCode postcheck(SNES,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
66: int main(int argc,char **argv)
67: {
68: SNES snes; /* nonlinear solver context */
69: Vec x,r; /* solution, residual vectors */
70: Mat J; /* Jacobian matrix */
71: AppCtx user; /* user-defined application context */
73: PetscInt i,its,N,hist_its[50];
74: PetscMPIInt size;
75: PetscReal bratu_lambda_max = 6.81,bratu_lambda_min = 0.,history[50];
76: MatFDColoring fdcoloring;
77: PetscBool matrix_free = PETSC_FALSE,flg,fd_coloring = PETSC_FALSE, use_convergence_test = PETSC_FALSE,pc = PETSC_FALSE;
78: KSP ksp;
79: PetscInt *testarray;
81: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
82: MPI_Comm_size(PETSC_COMM_WORLD,&size);
83: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"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(NULL,NULL,"-mx",&user.mx,NULL);
90: PetscOptionsGetInt(NULL,NULL,"-my",&user.my,NULL);
91: PetscOptionsGetReal(NULL,NULL,"-par",&user.param,NULL);
92: PetscOptionsGetBool(NULL,NULL,"-pc",&pc,NULL);
93: if (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lambda is out of range");
94: N = user.mx*user.my;
95: PetscOptionsGetBool(NULL,NULL,"-use_convergence_test",&use_convergence_test,NULL);
97: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
98: Create nonlinear solver context
99: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
101: SNESCreate(PETSC_COMM_WORLD,&snes);
103: if (pc) {
104: SNESSetType(snes,SNESNEWTONTR);
105: SNESNewtonTRSetPostCheck(snes, postcheck,NULL);
106: }
108: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109: Create vector data structures; set function evaluation routine
110: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
112: VecCreate(PETSC_COMM_WORLD,&x);
113: VecSetSizes(x,PETSC_DECIDE,N);
114: VecSetFromOptions(x);
115: VecDuplicate(x,&r);
117: /*
118: Set function evaluation routine and vector. Whenever the nonlinear
119: solver needs to evaluate the nonlinear function, it will call this
120: routine.
121: - Note that the final routine argument is the user-defined
122: context that provides application-specific data for the
123: function evaluation routine.
124: */
125: SNESSetFunction(snes,r,FormFunction,(void*)&user);
127: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: Create matrix data structure; set Jacobian evaluation routine
129: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
131: /*
132: Create matrix. Here we only approximately preallocate storage space
133: for the Jacobian. See the users manual for a discussion of better
134: techniques for preallocating matrix memory.
135: */
136: PetscOptionsGetBool(NULL,NULL,"-snes_mf",&matrix_free,NULL);
137: if (!matrix_free) {
138: PetscBool matrix_free_operator = PETSC_FALSE;
139: PetscOptionsGetBool(NULL,NULL,"-snes_mf_operator",&matrix_free_operator,NULL);
140: if (matrix_free_operator) matrix_free = PETSC_FALSE;
141: }
142: if (!matrix_free) {
143: MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,5,NULL,&J);
144: }
146: /*
147: This option will cause the Jacobian to be computed via finite differences
148: efficiently using a coloring of the columns of the matrix.
149: */
150: PetscOptionsGetBool(NULL,NULL,"-snes_fd_coloring",&fd_coloring,NULL);
151: if (matrix_free && fd_coloring) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_INCOMP,"Use only one of -snes_mf, -snes_fd_coloring options!\nYou can do -snes_mf_operator -snes_fd_coloring");
153: if (fd_coloring) {
154: ISColoring iscoloring;
155: MatColoring mc;
157: /*
158: This initializes the nonzero structure of the Jacobian. This is artificial
159: because clearly if we had a routine to compute the Jacobian we won't need
160: to use finite differences.
161: */
162: FormJacobian(snes,x,J,J,&user);
164: /*
165: Color the matrix, i.e. determine groups of columns that share no common
166: rows. These columns in the Jacobian can all be computed simultaneously.
167: */
168: MatColoringCreate(J,&mc);
169: MatColoringSetType(mc,MATCOLORINGSL);
170: MatColoringSetFromOptions(mc);
171: MatColoringApply(mc,&iscoloring);
172: MatColoringDestroy(&mc);
173: /*
174: Create the data structure that SNESComputeJacobianDefaultColor() uses
175: to compute the actual Jacobians via finite differences.
176: */
177: MatFDColoringCreate(J,iscoloring,&fdcoloring);
178: MatFDColoringSetFunction(fdcoloring,(PetscErrorCode (*)(void))FormFunction,&user);
179: MatFDColoringSetFromOptions(fdcoloring);
180: MatFDColoringSetUp(J,iscoloring,fdcoloring);
181: /*
182: Tell SNES to use the routine SNESComputeJacobianDefaultColor()
183: to compute Jacobians.
184: */
185: SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,fdcoloring);
186: ISColoringDestroy(&iscoloring);
187: }
188: /*
189: Set Jacobian matrix data structure and default Jacobian evaluation
190: routine. Whenever the nonlinear solver needs to compute the
191: Jacobian matrix, it will call this routine.
192: - Note that the final routine argument is the user-defined
193: context that provides application-specific data for the
194: Jacobian evaluation routine.
195: - The user can override with:
196: -snes_fd : default finite differencing approximation of Jacobian
197: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
198: (unless user explicitly sets preconditioner)
199: -snes_mf_operator : form preconditioning matrix as set by the user,
200: but use matrix-free approx for Jacobian-vector
201: products within Newton-Krylov method
202: */
203: else if (!matrix_free) {
204: SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);
205: }
207: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
208: Customize nonlinear solver; set runtime options
209: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
211: /*
212: Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
213: */
214: SNESSetFromOptions(snes);
216: /*
217: Set array that saves the function norms. This array is intended
218: when the user wants to save the convergence history for later use
219: rather than just to view the function norms via -snes_monitor.
220: */
221: SNESSetConvergenceHistory(snes,history,hist_its,50,PETSC_TRUE);
223: /*
224: Add a user provided convergence test; this is to test that SNESNEWTONTR properly calls the
225: user provided test before the specialized test. The convergence context is just an array to
226: test that it gets properly freed at the end
227: */
228: if (use_convergence_test) {
229: SNESGetKSP(snes,&ksp);
230: PetscMalloc1(5,&testarray);
231: KSPSetConvergenceTest(ksp,ConvergenceTest,testarray,ConvergenceDestroy);
232: }
234: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
235: Evaluate initial guess; then solve nonlinear system
236: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
237: /*
238: Note: The user should initialize the vector, x, with the initial guess
239: for the nonlinear solver prior to calling SNESSolve(). In particular,
240: to employ an initial guess of zero, the user should explicitly set
241: this vector to zero by calling VecSet().
242: */
243: FormInitialGuess(&user,x);
244: SNESSolve(snes,NULL,x);
245: SNESGetIterationNumber(snes,&its);
246: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
248: /*
249: Print the convergence history. This is intended just to demonstrate
250: use of the data attained via SNESSetConvergenceHistory().
251: */
252: PetscOptionsHasName(NULL,NULL,"-print_history",&flg);
253: if (flg) {
254: for (i=0; i<its+1; i++) {
255: PetscPrintf(PETSC_COMM_WORLD,"iteration %D: Linear iterations %D Function norm = %g\n",i,hist_its[i],(double)history[i]);
256: }
257: }
259: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
260: Free work space. All PETSc objects should be destroyed when they
261: are no longer needed.
262: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
264: if (!matrix_free) {
265: MatDestroy(&J);
266: }
267: if (fd_coloring) {
268: MatFDColoringDestroy(&fdcoloring);
269: }
270: VecDestroy(&x);
271: VecDestroy(&r);
272: SNESDestroy(&snes);
273: PetscFinalize();
274: return ierr;
275: }
276: /* ------------------------------------------------------------------- */
277: /*
278: FormInitialGuess - Forms initial approximation.
280: Input Parameters:
281: user - user-defined application context
282: X - vector
284: Output Parameter:
285: X - vector
286: */
287: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
288: {
289: PetscInt i,j,row,mx,my;
291: PetscReal lambda,temp1,temp,hx,hy;
292: PetscScalar *x;
294: mx = user->mx;
295: my = user->my;
296: lambda = user->param;
298: hx = 1.0 / (PetscReal)(mx-1);
299: hy = 1.0 / (PetscReal)(my-1);
301: /*
302: Get a pointer to vector data.
303: - For default PETSc vectors, VecGetArray() returns a pointer to
304: the data array. Otherwise, the routine is implementation dependent.
305: - You MUST call VecRestoreArray() when you no longer need access to
306: the array.
307: */
308: VecGetArray(X,&x);
309: temp1 = lambda/(lambda + 1.0);
310: for (j=0; j<my; j++) {
311: temp = (PetscReal)(PetscMin(j,my-j-1))*hy;
312: for (i=0; i<mx; i++) {
313: row = i + j*mx;
314: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
315: x[row] = 0.0;
316: continue;
317: }
318: x[row] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));
319: }
320: }
322: /*
323: Restore vector
324: */
325: VecRestoreArray(X,&x);
326: return 0;
327: }
328: /* ------------------------------------------------------------------- */
329: /*
330: FormFunction - Evaluates nonlinear function, F(x).
332: Input Parameters:
333: . snes - the SNES context
334: . X - input vector
335: . ptr - optional user-defined context, as set by SNESSetFunction()
337: Output Parameter:
338: . F - function vector
339: */
340: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
341: {
342: AppCtx *user = (AppCtx*)ptr;
343: PetscInt i,j,row,mx,my;
344: PetscErrorCode ierr;
345: PetscReal two = 2.0,one = 1.0,lambda,hx,hy,hxdhy,hydhx;
346: PetscScalar ut,ub,ul,ur,u,uxx,uyy,sc,*f;
347: const PetscScalar *x;
349: mx = user->mx;
350: my = user->my;
351: lambda = user->param;
352: hx = one / (PetscReal)(mx-1);
353: hy = one / (PetscReal)(my-1);
354: sc = hx*hy;
355: hxdhy = hx/hy;
356: hydhx = hy/hx;
358: /*
359: Get pointers to vector data
360: */
361: VecGetArrayRead(X,&x);
362: VecGetArray(F,&f);
364: /*
365: Compute function
366: */
367: for (j=0; j<my; j++) {
368: for (i=0; i<mx; i++) {
369: row = i + j*mx;
370: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
371: f[row] = x[row];
372: continue;
373: }
374: u = x[row];
375: ub = x[row - mx];
376: ul = x[row - 1];
377: ut = x[row + mx];
378: ur = x[row + 1];
379: uxx = (-ur + two*u - ul)*hydhx;
380: uyy = (-ut + two*u - ub)*hxdhy;
381: f[row] = uxx + uyy - sc*lambda*PetscExpScalar(u);
382: }
383: }
385: /*
386: Restore vectors
387: */
388: VecRestoreArrayRead(X,&x);
389: VecRestoreArray(F,&f);
390: return 0;
391: }
392: /* ------------------------------------------------------------------- */
393: /*
394: FormJacobian - Evaluates Jacobian matrix.
396: Input Parameters:
397: . snes - the SNES context
398: . x - input vector
399: . ptr - optional user-defined context, as set by SNESSetJacobian()
401: Output Parameters:
402: . A - Jacobian matrix
403: . B - optionally different preconditioning matrix
404: . flag - flag indicating matrix structure
405: */
406: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat jac,void *ptr)
407: {
408: AppCtx *user = (AppCtx*)ptr; /* user-defined applicatin context */
409: PetscInt i,j,row,mx,my,col[5];
410: PetscErrorCode ierr;
411: PetscScalar two = 2.0,one = 1.0,lambda,v[5],sc;
412: const PetscScalar *x;
413: PetscReal hx,hy,hxdhy,hydhx;
415: mx = user->mx;
416: my = user->my;
417: lambda = user->param;
418: hx = 1.0 / (PetscReal)(mx-1);
419: hy = 1.0 / (PetscReal)(my-1);
420: sc = hx*hy;
421: hxdhy = hx/hy;
422: hydhx = hy/hx;
424: /*
425: Get pointer to vector data
426: */
427: VecGetArrayRead(X,&x);
429: /*
430: Compute entries of the Jacobian
431: */
432: for (j=0; j<my; j++) {
433: for (i=0; i<mx; i++) {
434: row = i + j*mx;
435: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
436: MatSetValues(jac,1,&row,1,&row,&one,INSERT_VALUES);
437: continue;
438: }
439: v[0] = -hxdhy; col[0] = row - mx;
440: v[1] = -hydhx; col[1] = row - 1;
441: v[2] = two*(hydhx + hxdhy) - sc*lambda*PetscExpScalar(x[row]); col[2] = row;
442: v[3] = -hydhx; col[3] = row + 1;
443: v[4] = -hxdhy; col[4] = row + mx;
444: MatSetValues(jac,1,&row,5,col,v,INSERT_VALUES);
445: }
446: }
448: /*
449: Restore vector
450: */
451: VecRestoreArrayRead(X,&x);
453: /*
454: Assemble matrix
455: */
456: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
457: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
459: if (jac != J) {
460: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
461: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
462: }
464: return 0;
465: }
467: PetscErrorCode ConvergenceTest(KSP ksp,PetscInt it,PetscReal nrm,KSPConvergedReason *reason,void *ctx)
468: {
472: *reason = KSP_CONVERGED_ITERATING;
473: if (it > 1) {
474: *reason = KSP_CONVERGED_ITS;
475: PetscInfo(NULL,"User provided convergence test returning after 2 iterations\n");
476: }
477: return(0);
478: }
480: PetscErrorCode ConvergenceDestroy(void* ctx)
481: {
485: PetscInfo(NULL,"User provided convergence destroy called\n");
486: PetscFree(ctx);
487: return(0);
488: }
490: PetscErrorCode postcheck(SNES snes,Vec x,Vec y,Vec w,PetscBool *changed_y,PetscBool *changed_w,void *ctx)
491: {
493: PetscReal norm;
494: Vec tmp;
497: VecDuplicate(x,&tmp);
498: VecWAXPY(tmp,-1.0,x,w);
499: VecNorm(tmp,NORM_2,&norm);
500: VecDestroy(&tmp);
501: PetscPrintf(PETSC_COMM_WORLD,"Norm of search step %g\n",(double)norm);
502: return(0);
503: }
505: /*TEST
507: build:
508: requires: !single
510: test:
511: args: -ksp_gmres_cgs_refinement_type refine_always
513: test:
514: suffix: 2
515: args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always
517: test:
518: suffix: 2a
519: filter: grep -i KSPConvergedDefault > /dev/null && echo "Found KSPConvergedDefault"
520: args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always -info
521: requires: defined(PETSC_USE_INFO)
523: test:
524: suffix: 2b
525: filter: grep -i "User provided convergence test" > /dev/null && echo "Found User provided convergence test"
526: args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always -use_convergence_test -info
527: requires: defined(PETSC_USE_INFO)
529: test:
530: suffix: 3
531: args: -snes_monitor_short -mat_coloring_type sl -snes_fd_coloring -mx 8 -my 11 -ksp_gmres_cgs_refinement_type refine_always
533: test:
534: suffix: 4
535: args: -pc -par 6.807 -snes_monitor -snes_converged_reason
537: TEST*/