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 */
72: PetscInt i,its,N,hist_its[50];
73: PetscMPIInt size;
74: PetscReal bratu_lambda_max = 6.81,bratu_lambda_min = 0.,history[50];
75: MatFDColoring fdcoloring;
76: PetscBool matrix_free = PETSC_FALSE,flg,fd_coloring = PETSC_FALSE, use_convergence_test = PETSC_FALSE,pc = PETSC_FALSE;
77: KSP ksp;
78: PetscInt *testarray;
80: PetscInitialize(&argc,&argv,(char*)0,help);
81: MPI_Comm_size(PETSC_COMM_WORLD,&size);
84: /*
85: Initialize problem parameters
86: */
87: user.mx = 4; user.my = 4; user.param = 6.0;
88: PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,NULL);
89: PetscOptionsGetInt(NULL,NULL,"-my",&user.my,NULL);
90: PetscOptionsGetReal(NULL,NULL,"-par",&user.param,NULL);
91: PetscOptionsGetBool(NULL,NULL,"-pc",&pc,NULL);
93: N = user.mx*user.my;
94: PetscOptionsGetBool(NULL,NULL,"-use_convergence_test",&use_convergence_test,NULL);
96: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97: Create nonlinear solver context
98: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
100: SNESCreate(PETSC_COMM_WORLD,&snes);
102: if (pc) {
103: SNESSetType(snes,SNESNEWTONTR);
104: SNESNewtonTRSetPostCheck(snes, postcheck,NULL);
105: }
107: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108: Create vector data structures; set function evaluation routine
109: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
111: VecCreate(PETSC_COMM_WORLD,&x);
112: VecSetSizes(x,PETSC_DECIDE,N);
113: VecSetFromOptions(x);
114: VecDuplicate(x,&r);
116: /*
117: Set function evaluation routine and vector. Whenever the nonlinear
118: solver needs to evaluate the nonlinear function, it will call this
119: routine.
120: - Note that the final routine argument is the user-defined
121: context that provides application-specific data for the
122: function evaluation routine.
123: */
124: SNESSetFunction(snes,r,FormFunction,(void*)&user);
126: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127: Create matrix data structure; set Jacobian evaluation routine
128: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130: /*
131: Create matrix. Here we only approximately preallocate storage space
132: for the Jacobian. See the users manual for a discussion of better
133: techniques for preallocating matrix memory.
134: */
135: PetscOptionsGetBool(NULL,NULL,"-snes_mf",&matrix_free,NULL);
136: if (!matrix_free) {
137: PetscBool matrix_free_operator = PETSC_FALSE;
138: PetscOptionsGetBool(NULL,NULL,"-snes_mf_operator",&matrix_free_operator,NULL);
139: if (matrix_free_operator) matrix_free = PETSC_FALSE;
140: }
141: if (!matrix_free) {
142: MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,5,NULL,&J);
143: }
145: /*
146: This option will cause the Jacobian to be computed via finite differences
147: efficiently using a coloring of the columns of the matrix.
148: */
149: PetscOptionsGetBool(NULL,NULL,"-snes_fd_coloring",&fd_coloring,NULL);
152: if (fd_coloring) {
153: ISColoring iscoloring;
154: MatColoring mc;
156: /*
157: This initializes the nonzero structure of the Jacobian. This is artificial
158: because clearly if we had a routine to compute the Jacobian we won't need
159: to use finite differences.
160: */
161: FormJacobian(snes,x,J,J,&user);
163: /*
164: Color the matrix, i.e. determine groups of columns that share no common
165: rows. These columns in the Jacobian can all be computed simultaneously.
166: */
167: MatColoringCreate(J,&mc);
168: MatColoringSetType(mc,MATCOLORINGSL);
169: MatColoringSetFromOptions(mc);
170: MatColoringApply(mc,&iscoloring);
171: MatColoringDestroy(&mc);
172: /*
173: Create the data structure that SNESComputeJacobianDefaultColor() uses
174: to compute the actual Jacobians via finite differences.
175: */
176: MatFDColoringCreate(J,iscoloring,&fdcoloring);
177: MatFDColoringSetFunction(fdcoloring,(PetscErrorCode (*)(void))FormFunction,&user);
178: MatFDColoringSetFromOptions(fdcoloring);
179: MatFDColoringSetUp(J,iscoloring,fdcoloring);
180: /*
181: Tell SNES to use the routine SNESComputeJacobianDefaultColor()
182: to compute Jacobians.
183: */
184: SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,fdcoloring);
185: ISColoringDestroy(&iscoloring);
186: }
187: /*
188: Set Jacobian matrix data structure and default Jacobian evaluation
189: routine. Whenever the nonlinear solver needs to compute the
190: Jacobian matrix, it will call this routine.
191: - Note that the final routine argument is the user-defined
192: context that provides application-specific data for the
193: Jacobian evaluation routine.
194: - The user can override with:
195: -snes_fd : default finite differencing approximation of Jacobian
196: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
197: (unless user explicitly sets preconditioner)
198: -snes_mf_operator : form preconditioning matrix as set by the user,
199: but use matrix-free approx for Jacobian-vector
200: products within Newton-Krylov method
201: */
202: else if (!matrix_free) {
203: SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);
204: }
206: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
207: Customize nonlinear solver; set runtime options
208: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
210: /*
211: Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
212: */
213: SNESSetFromOptions(snes);
215: /*
216: Set array that saves the function norms. This array is intended
217: when the user wants to save the convergence history for later use
218: rather than just to view the function norms via -snes_monitor.
219: */
220: SNESSetConvergenceHistory(snes,history,hist_its,50,PETSC_TRUE);
222: /*
223: Add a user provided convergence test; this is to test that SNESNEWTONTR properly calls the
224: user provided test before the specialized test. The convergence context is just an array to
225: test that it gets properly freed at the end
226: */
227: if (use_convergence_test) {
228: SNESGetKSP(snes,&ksp);
229: PetscMalloc1(5,&testarray);
230: KSPSetConvergenceTest(ksp,ConvergenceTest,testarray,ConvergenceDestroy);
231: }
233: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
234: Evaluate initial guess; then solve nonlinear system
235: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
236: /*
237: Note: The user should initialize the vector, x, with the initial guess
238: for the nonlinear solver prior to calling SNESSolve(). In particular,
239: to employ an initial guess of zero, the user should explicitly set
240: this vector to zero by calling VecSet().
241: */
242: FormInitialGuess(&user,x);
243: SNESSolve(snes,NULL,x);
244: SNESGetIterationNumber(snes,&its);
245: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
247: /*
248: Print the convergence history. This is intended just to demonstrate
249: use of the data attained via SNESSetConvergenceHistory().
250: */
251: PetscOptionsHasName(NULL,NULL,"-print_history",&flg);
252: if (flg) {
253: for (i=0; i<its+1; i++) {
254: PetscPrintf(PETSC_COMM_WORLD,"iteration %D: Linear iterations %D Function norm = %g\n",i,hist_its[i],(double)history[i]);
255: }
256: }
258: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
259: Free work space. All PETSc objects should be destroyed when they
260: are no longer needed.
261: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
263: if (!matrix_free) {
264: MatDestroy(&J);
265: }
266: if (fd_coloring) {
267: MatFDColoringDestroy(&fdcoloring);
268: }
269: VecDestroy(&x);
270: VecDestroy(&r);
271: SNESDestroy(&snes);
272: PetscFinalize();
273: return 0;
274: }
275: /* ------------------------------------------------------------------- */
276: /*
277: FormInitialGuess - Forms initial approximation.
279: Input Parameters:
280: user - user-defined application context
281: X - vector
283: Output Parameter:
284: X - vector
285: */
286: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
287: {
288: PetscInt i,j,row,mx,my;
289: PetscReal lambda,temp1,temp,hx,hy;
290: PetscScalar *x;
292: mx = user->mx;
293: my = user->my;
294: lambda = user->param;
296: hx = 1.0 / (PetscReal)(mx-1);
297: hy = 1.0 / (PetscReal)(my-1);
299: /*
300: Get a pointer to vector data.
301: - For default PETSc vectors, VecGetArray() returns a pointer to
302: the data array. Otherwise, the routine is implementation dependent.
303: - You MUST call VecRestoreArray() when you no longer need access to
304: the array.
305: */
306: VecGetArray(X,&x);
307: temp1 = lambda/(lambda + 1.0);
308: for (j=0; j<my; j++) {
309: temp = (PetscReal)(PetscMin(j,my-j-1))*hy;
310: for (i=0; i<mx; i++) {
311: row = i + j*mx;
312: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
313: x[row] = 0.0;
314: continue;
315: }
316: x[row] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));
317: }
318: }
320: /*
321: Restore vector
322: */
323: VecRestoreArray(X,&x);
324: return 0;
325: }
326: /* ------------------------------------------------------------------- */
327: /*
328: FormFunction - Evaluates nonlinear function, F(x).
330: Input Parameters:
331: . snes - the SNES context
332: . X - input vector
333: . ptr - optional user-defined context, as set by SNESSetFunction()
335: Output Parameter:
336: . F - function vector
337: */
338: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
339: {
340: AppCtx *user = (AppCtx*)ptr;
341: PetscInt i,j,row,mx,my;
342: PetscReal two = 2.0,one = 1.0,lambda,hx,hy,hxdhy,hydhx;
343: PetscScalar ut,ub,ul,ur,u,uxx,uyy,sc,*f;
344: const PetscScalar *x;
346: mx = user->mx;
347: my = user->my;
348: lambda = user->param;
349: hx = one / (PetscReal)(mx-1);
350: hy = one / (PetscReal)(my-1);
351: sc = hx*hy;
352: hxdhy = hx/hy;
353: hydhx = hy/hx;
355: /*
356: Get pointers to vector data
357: */
358: VecGetArrayRead(X,&x);
359: VecGetArray(F,&f);
361: /*
362: Compute function
363: */
364: for (j=0; j<my; j++) {
365: for (i=0; i<mx; i++) {
366: row = i + j*mx;
367: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
368: f[row] = x[row];
369: continue;
370: }
371: u = x[row];
372: ub = x[row - mx];
373: ul = x[row - 1];
374: ut = x[row + mx];
375: ur = x[row + 1];
376: uxx = (-ur + two*u - ul)*hydhx;
377: uyy = (-ut + two*u - ub)*hxdhy;
378: f[row] = uxx + uyy - sc*lambda*PetscExpScalar(u);
379: }
380: }
382: /*
383: Restore vectors
384: */
385: VecRestoreArrayRead(X,&x);
386: VecRestoreArray(F,&f);
387: return 0;
388: }
389: /* ------------------------------------------------------------------- */
390: /*
391: FormJacobian - Evaluates Jacobian matrix.
393: Input Parameters:
394: . snes - the SNES context
395: . x - input vector
396: . ptr - optional user-defined context, as set by SNESSetJacobian()
398: Output Parameters:
399: . A - Jacobian matrix
400: . B - optionally different preconditioning matrix
401: . flag - flag indicating matrix structure
402: */
403: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat jac,void *ptr)
404: {
405: AppCtx *user = (AppCtx*)ptr; /* user-defined applicatin context */
406: PetscInt i,j,row,mx,my,col[5];
407: PetscScalar two = 2.0,one = 1.0,lambda,v[5],sc;
408: const PetscScalar *x;
409: PetscReal hx,hy,hxdhy,hydhx;
411: mx = user->mx;
412: my = user->my;
413: lambda = user->param;
414: hx = 1.0 / (PetscReal)(mx-1);
415: hy = 1.0 / (PetscReal)(my-1);
416: sc = hx*hy;
417: hxdhy = hx/hy;
418: hydhx = hy/hx;
420: /*
421: Get pointer to vector data
422: */
423: VecGetArrayRead(X,&x);
425: /*
426: Compute entries of the Jacobian
427: */
428: for (j=0; j<my; j++) {
429: for (i=0; i<mx; i++) {
430: row = i + j*mx;
431: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
432: MatSetValues(jac,1,&row,1,&row,&one,INSERT_VALUES);
433: continue;
434: }
435: v[0] = -hxdhy; col[0] = row - mx;
436: v[1] = -hydhx; col[1] = row - 1;
437: v[2] = two*(hydhx + hxdhy) - sc*lambda*PetscExpScalar(x[row]); col[2] = row;
438: v[3] = -hydhx; col[3] = row + 1;
439: v[4] = -hxdhy; col[4] = row + mx;
440: MatSetValues(jac,1,&row,5,col,v,INSERT_VALUES);
441: }
442: }
444: /*
445: Restore vector
446: */
447: VecRestoreArrayRead(X,&x);
449: /*
450: Assemble matrix
451: */
452: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
453: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
455: if (jac != J) {
456: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
457: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
458: }
460: return 0;
461: }
463: PetscErrorCode ConvergenceTest(KSP ksp,PetscInt it,PetscReal nrm,KSPConvergedReason *reason,void *ctx)
464: {
465: *reason = KSP_CONVERGED_ITERATING;
466: if (it > 1) {
467: *reason = KSP_CONVERGED_ITS;
468: PetscInfo(NULL,"User provided convergence test returning after 2 iterations\n");
469: }
470: return 0;
471: }
473: PetscErrorCode ConvergenceDestroy(void* ctx)
474: {
475: PetscInfo(NULL,"User provided convergence destroy called\n");
476: PetscFree(ctx);
477: return 0;
478: }
480: PetscErrorCode postcheck(SNES snes,Vec x,Vec y,Vec w,PetscBool *changed_y,PetscBool *changed_w,void *ctx)
481: {
482: PetscReal norm;
483: Vec tmp;
485: VecDuplicate(x,&tmp);
486: VecWAXPY(tmp,-1.0,x,w);
487: VecNorm(tmp,NORM_2,&norm);
488: VecDestroy(&tmp);
489: PetscPrintf(PETSC_COMM_WORLD,"Norm of search step %g\n",(double)norm);
490: return 0;
491: }
493: /*TEST
495: build:
496: requires: !single
498: test:
499: args: -ksp_gmres_cgs_refinement_type refine_always
501: test:
502: suffix: 2
503: args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always
505: test:
506: suffix: 2a
507: filter: grep -i KSPConvergedDefault > /dev/null && echo "Found KSPConvergedDefault"
508: args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always -info
509: requires: defined(PETSC_USE_INFO)
511: test:
512: suffix: 2b
513: filter: grep -i "User provided convergence test" > /dev/null && echo "Found User provided convergence test"
514: args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always -use_convergence_test -info
515: requires: defined(PETSC_USE_INFO)
517: test:
518: suffix: 3
519: args: -snes_monitor_short -mat_coloring_type sl -snes_fd_coloring -mx 8 -my 11 -ksp_gmres_cgs_refinement_type refine_always
521: test:
522: suffix: 4
523: args: -pc -par 6.807 -snes_monitor -snes_converged_reason
525: TEST*/