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: /* ------------------------------------------------------------------------
11: Solid Fuel Ignition (SFI) problem. This problem is modeled by
12: the partial differential equation
14: -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
16: with boundary conditions
18: u = 0 for x = 0, x = 1, y = 0, y = 1.
20: A finite difference approximation with the usual 5-point stencil
21: is used to discretize the boundary value problem to obtain a nonlinear
22: system of equations.
24: The parallel version of this code is snes/tutorials/ex5.c
26: ------------------------------------------------------------------------- */
28: /*
29: Include "petscsnes.h" so that we can use SNES solvers. Note that
30: this file automatically includes:
31: petscsys.h - base PETSc routines petscvec.h - vectors
32: petscmat.h - matrices
33: petscis.h - index sets petscksp.h - Krylov subspace methods
34: petscviewer.h - viewers petscpc.h - preconditioners
35: petscksp.h - linear solvers
36: */
38: #include <petscsnes.h>
40: /*
41: User-defined application context - contains data needed by the
42: application-provided call-back routines, FormJacobian() and
43: FormFunction().
44: */
45: typedef struct {
46: PetscReal param; /* test problem parameter */
47: PetscInt mx; /* Discretization in x-direction */
48: PetscInt my; /* Discretization in y-direction */
49: } AppCtx;
51: /*
52: User-defined routines
53: */
54: extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
55: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
56: extern PetscErrorCode FormInitialGuess(AppCtx*,Vec);
57: extern PetscErrorCode ConvergenceTest(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*);
58: extern PetscErrorCode ConvergenceDestroy(void*);
59: extern PetscErrorCode postcheck(SNES,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
61: int main(int argc,char **argv)
62: {
63: SNES snes; /* nonlinear solver context */
64: Vec x,r; /* solution, residual vectors */
65: Mat J; /* Jacobian matrix */
66: AppCtx user; /* user-defined application context */
67: PetscInt i,its,N,hist_its[50];
68: PetscMPIInt size;
69: PetscReal bratu_lambda_max = 6.81,bratu_lambda_min = 0.,history[50];
70: MatFDColoring fdcoloring;
71: PetscBool matrix_free = PETSC_FALSE,flg,fd_coloring = PETSC_FALSE, use_convergence_test = PETSC_FALSE,pc = PETSC_FALSE;
72: KSP ksp;
73: PetscInt *testarray;
75: PetscInitialize(&argc,&argv,(char*)0,help);
76: MPI_Comm_size(PETSC_COMM_WORLD,&size);
79: /*
80: Initialize problem parameters
81: */
82: user.mx = 4; user.my = 4; user.param = 6.0;
83: PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,NULL);
84: PetscOptionsGetInt(NULL,NULL,"-my",&user.my,NULL);
85: PetscOptionsGetReal(NULL,NULL,"-par",&user.param,NULL);
86: PetscOptionsGetBool(NULL,NULL,"-pc",&pc,NULL);
88: N = user.mx*user.my;
89: PetscOptionsGetBool(NULL,NULL,"-use_convergence_test",&use_convergence_test,NULL);
91: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
92: Create nonlinear solver context
93: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
95: SNESCreate(PETSC_COMM_WORLD,&snes);
97: if (pc) {
98: SNESSetType(snes,SNESNEWTONTR);
99: SNESNewtonTRSetPostCheck(snes, postcheck,NULL);
100: }
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: PetscOptionsGetBool(NULL,NULL,"-snes_mf",&matrix_free,NULL);
131: if (!matrix_free) {
132: PetscBool matrix_free_operator = PETSC_FALSE;
133: PetscOptionsGetBool(NULL,NULL,"-snes_mf_operator",&matrix_free_operator,NULL);
134: if (matrix_free_operator) matrix_free = PETSC_FALSE;
135: }
136: if (!matrix_free) {
137: MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,5,NULL,&J);
138: }
140: /*
141: This option will cause the Jacobian to be computed via finite differences
142: efficiently using a coloring of the columns of the matrix.
143: */
144: PetscOptionsGetBool(NULL,NULL,"-snes_fd_coloring",&fd_coloring,NULL);
147: if (fd_coloring) {
148: ISColoring iscoloring;
149: MatColoring mc;
151: /*
152: This initializes the nonzero structure of the Jacobian. This is artificial
153: because clearly if we had a routine to compute the Jacobian we won't need
154: to use finite differences.
155: */
156: FormJacobian(snes,x,J,J,&user);
158: /*
159: Color the matrix, i.e. determine groups of columns that share no common
160: rows. These columns in the Jacobian can all be computed simultaneously.
161: */
162: MatColoringCreate(J,&mc);
163: MatColoringSetType(mc,MATCOLORINGSL);
164: MatColoringSetFromOptions(mc);
165: MatColoringApply(mc,&iscoloring);
166: MatColoringDestroy(&mc);
167: /*
168: Create the data structure that SNESComputeJacobianDefaultColor() uses
169: to compute the actual Jacobians via finite differences.
170: */
171: MatFDColoringCreate(J,iscoloring,&fdcoloring);
172: MatFDColoringSetFunction(fdcoloring,(PetscErrorCode (*)(void))FormFunction,&user);
173: MatFDColoringSetFromOptions(fdcoloring);
174: MatFDColoringSetUp(J,iscoloring,fdcoloring);
175: /*
176: Tell SNES to use the routine SNESComputeJacobianDefaultColor()
177: to compute Jacobians.
178: */
179: SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,fdcoloring);
180: ISColoringDestroy(&iscoloring);
181: }
182: /*
183: Set Jacobian matrix data structure and default Jacobian evaluation
184: routine. Whenever the nonlinear solver needs to compute the
185: Jacobian matrix, it will call this routine.
186: - Note that the final routine argument is the user-defined
187: context that provides application-specific data for the
188: Jacobian evaluation routine.
189: - The user can override with:
190: -snes_fd : default finite differencing approximation of Jacobian
191: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
192: (unless user explicitly sets preconditioner)
193: -snes_mf_operator : form preconditioning matrix as set by the user,
194: but use matrix-free approx for Jacobian-vector
195: products within Newton-Krylov method
196: */
197: else if (!matrix_free) {
198: SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);
199: }
201: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202: Customize nonlinear solver; set runtime options
203: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
205: /*
206: Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
207: */
208: SNESSetFromOptions(snes);
210: /*
211: Set array that saves the function norms. This array is intended
212: when the user wants to save the convergence history for later use
213: rather than just to view the function norms via -snes_monitor.
214: */
215: SNESSetConvergenceHistory(snes,history,hist_its,50,PETSC_TRUE);
217: /*
218: Add a user provided convergence test; this is to test that SNESNEWTONTR properly calls the
219: user provided test before the specialized test. The convergence context is just an array to
220: test that it gets properly freed at the end
221: */
222: if (use_convergence_test) {
223: SNESGetKSP(snes,&ksp);
224: PetscMalloc1(5,&testarray);
225: KSPSetConvergenceTest(ksp,ConvergenceTest,testarray,ConvergenceDestroy);
226: }
228: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
229: Evaluate initial guess; then solve nonlinear system
230: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
231: /*
232: Note: The user should initialize the vector, x, with the initial guess
233: for the nonlinear solver prior to calling SNESSolve(). In particular,
234: to employ an initial guess of zero, the user should explicitly set
235: this vector to zero by calling VecSet().
236: */
237: FormInitialGuess(&user,x);
238: SNESSolve(snes,NULL,x);
239: SNESGetIterationNumber(snes,&its);
240: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
242: /*
243: Print the convergence history. This is intended just to demonstrate
244: use of the data attained via SNESSetConvergenceHistory().
245: */
246: PetscOptionsHasName(NULL,NULL,"-print_history",&flg);
247: if (flg) {
248: for (i=0; i<its+1; i++) {
249: PetscPrintf(PETSC_COMM_WORLD,"iteration %D: Linear iterations %D Function norm = %g\n",i,hist_its[i],(double)history[i]);
250: }
251: }
253: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
254: Free work space. All PETSc objects should be destroyed when they
255: are no longer needed.
256: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
258: if (!matrix_free) {
259: MatDestroy(&J);
260: }
261: if (fd_coloring) {
262: MatFDColoringDestroy(&fdcoloring);
263: }
264: VecDestroy(&x);
265: VecDestroy(&r);
266: SNESDestroy(&snes);
267: PetscFinalize();
268: return 0;
269: }
270: /* ------------------------------------------------------------------- */
271: /*
272: FormInitialGuess - Forms initial approximation.
274: Input Parameters:
275: user - user-defined application context
276: X - vector
278: Output Parameter:
279: X - vector
280: */
281: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
282: {
283: PetscInt i,j,row,mx,my;
284: PetscReal lambda,temp1,temp,hx,hy;
285: PetscScalar *x;
287: mx = user->mx;
288: my = user->my;
289: lambda = user->param;
291: hx = 1.0 / (PetscReal)(mx-1);
292: hy = 1.0 / (PetscReal)(my-1);
294: /*
295: Get a pointer to vector data.
296: - For default PETSc vectors, VecGetArray() returns a pointer to
297: the data array. Otherwise, the routine is implementation dependent.
298: - You MUST call VecRestoreArray() when you no longer need access to
299: the array.
300: */
301: VecGetArray(X,&x);
302: temp1 = lambda/(lambda + 1.0);
303: for (j=0; j<my; j++) {
304: temp = (PetscReal)(PetscMin(j,my-j-1))*hy;
305: for (i=0; i<mx; i++) {
306: row = i + j*mx;
307: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
308: x[row] = 0.0;
309: continue;
310: }
311: x[row] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));
312: }
313: }
315: /*
316: Restore vector
317: */
318: VecRestoreArray(X,&x);
319: return 0;
320: }
321: /* ------------------------------------------------------------------- */
322: /*
323: FormFunction - Evaluates nonlinear function, F(x).
325: Input Parameters:
326: . snes - the SNES context
327: . X - input vector
328: . ptr - optional user-defined context, as set by SNESSetFunction()
330: Output Parameter:
331: . F - function vector
332: */
333: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
334: {
335: AppCtx *user = (AppCtx*)ptr;
336: PetscInt i,j,row,mx,my;
337: PetscReal two = 2.0,one = 1.0,lambda,hx,hy,hxdhy,hydhx;
338: PetscScalar ut,ub,ul,ur,u,uxx,uyy,sc,*f;
339: const PetscScalar *x;
341: mx = user->mx;
342: my = user->my;
343: lambda = user->param;
344: hx = one / (PetscReal)(mx-1);
345: hy = one / (PetscReal)(my-1);
346: sc = hx*hy;
347: hxdhy = hx/hy;
348: hydhx = hy/hx;
350: /*
351: Get pointers to vector data
352: */
353: VecGetArrayRead(X,&x);
354: VecGetArray(F,&f);
356: /*
357: Compute function
358: */
359: for (j=0; j<my; j++) {
360: for (i=0; i<mx; i++) {
361: row = i + j*mx;
362: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
363: f[row] = x[row];
364: continue;
365: }
366: u = x[row];
367: ub = x[row - mx];
368: ul = x[row - 1];
369: ut = x[row + mx];
370: ur = x[row + 1];
371: uxx = (-ur + two*u - ul)*hydhx;
372: uyy = (-ut + two*u - ub)*hxdhy;
373: f[row] = uxx + uyy - sc*lambda*PetscExpScalar(u);
374: }
375: }
377: /*
378: Restore vectors
379: */
380: VecRestoreArrayRead(X,&x);
381: VecRestoreArray(F,&f);
382: return 0;
383: }
384: /* ------------------------------------------------------------------- */
385: /*
386: FormJacobian - Evaluates Jacobian matrix.
388: Input Parameters:
389: . snes - the SNES context
390: . x - input vector
391: . ptr - optional user-defined context, as set by SNESSetJacobian()
393: Output Parameters:
394: . A - Jacobian matrix
395: . B - optionally different preconditioning matrix
396: . flag - flag indicating matrix structure
397: */
398: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat jac,void *ptr)
399: {
400: AppCtx *user = (AppCtx*)ptr; /* user-defined applicatin context */
401: PetscInt i,j,row,mx,my,col[5];
402: PetscScalar two = 2.0,one = 1.0,lambda,v[5],sc;
403: const PetscScalar *x;
404: PetscReal hx,hy,hxdhy,hydhx;
406: mx = user->mx;
407: my = user->my;
408: lambda = user->param;
409: hx = 1.0 / (PetscReal)(mx-1);
410: hy = 1.0 / (PetscReal)(my-1);
411: sc = hx*hy;
412: hxdhy = hx/hy;
413: hydhx = hy/hx;
415: /*
416: Get pointer to vector data
417: */
418: VecGetArrayRead(X,&x);
420: /*
421: Compute entries of the Jacobian
422: */
423: for (j=0; j<my; j++) {
424: for (i=0; i<mx; i++) {
425: row = i + j*mx;
426: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
427: MatSetValues(jac,1,&row,1,&row,&one,INSERT_VALUES);
428: continue;
429: }
430: v[0] = -hxdhy; col[0] = row - mx;
431: v[1] = -hydhx; col[1] = row - 1;
432: v[2] = two*(hydhx + hxdhy) - sc*lambda*PetscExpScalar(x[row]); col[2] = row;
433: v[3] = -hydhx; col[3] = row + 1;
434: v[4] = -hxdhy; col[4] = row + mx;
435: MatSetValues(jac,1,&row,5,col,v,INSERT_VALUES);
436: }
437: }
439: /*
440: Restore vector
441: */
442: VecRestoreArrayRead(X,&x);
444: /*
445: Assemble matrix
446: */
447: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
448: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
450: if (jac != J) {
451: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
452: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
453: }
455: return 0;
456: }
458: PetscErrorCode ConvergenceTest(KSP ksp,PetscInt it,PetscReal nrm,KSPConvergedReason *reason,void *ctx)
459: {
460: *reason = KSP_CONVERGED_ITERATING;
461: if (it > 1) {
462: *reason = KSP_CONVERGED_ITS;
463: PetscInfo(NULL,"User provided convergence test returning after 2 iterations\n");
464: }
465: return 0;
466: }
468: PetscErrorCode ConvergenceDestroy(void* ctx)
469: {
470: PetscInfo(NULL,"User provided convergence destroy called\n");
471: PetscFree(ctx);
472: return 0;
473: }
475: PetscErrorCode postcheck(SNES snes,Vec x,Vec y,Vec w,PetscBool *changed_y,PetscBool *changed_w,void *ctx)
476: {
477: PetscReal norm;
478: Vec tmp;
480: VecDuplicate(x,&tmp);
481: VecWAXPY(tmp,-1.0,x,w);
482: VecNorm(tmp,NORM_2,&norm);
483: VecDestroy(&tmp);
484: PetscPrintf(PETSC_COMM_WORLD,"Norm of search step %g\n",(double)norm);
485: return 0;
486: }
488: /*TEST
490: build:
491: requires: !single
493: test:
494: args: -ksp_gmres_cgs_refinement_type refine_always
496: test:
497: suffix: 2
498: args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always
500: test:
501: suffix: 2a
502: filter: grep -i KSPConvergedDefault > /dev/null && echo "Found KSPConvergedDefault"
503: args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always -info
504: requires: defined(PETSC_USE_INFO)
506: test:
507: suffix: 2b
508: filter: grep -i "User provided convergence test" > /dev/null && echo "Found User provided convergence test"
509: args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always -use_convergence_test -info
510: requires: defined(PETSC_USE_INFO)
512: test:
513: suffix: 3
514: args: -snes_monitor_short -mat_coloring_type sl -snes_fd_coloring -mx 8 -my 11 -ksp_gmres_cgs_refinement_type refine_always
516: test:
517: suffix: 4
518: args: -pc -par 6.807 -snes_monitor -snes_converged_reason
520: TEST*/