Actual source code: ex1.c
petsc-3.11.3 2019-06-26
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*/
16: /* ------------------------------------------------------------------------
18: Solid Fuel Ignition (SFI) problem. This problem is modeled by
19: the partial differential equation
21: -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
23: with boundary conditions
25: u = 0 for x = 0, x = 1, y = 0, y = 1.
27: A finite difference approximation with the usual 5-point stencil
28: is used to discretize the boundary value problem to obtain a nonlinear
29: system of equations.
31: The parallel version of this code is snes/examples/tutorials/ex5.c
33: ------------------------------------------------------------------------- */
35: /*
36: Include "petscsnes.h" so that we can use SNES solvers. Note that
37: this file automatically includes:
38: petscsys.h - base PETSc routines petscvec.h - vectors
39: petscmat.h - matrices
40: petscis.h - index sets petscksp.h - Krylov subspace methods
41: petscviewer.h - viewers petscpc.h - preconditioners
42: petscksp.h - linear solvers
43: */
45: #include <petscsnes.h>
47: /*
48: User-defined application context - contains data needed by the
49: application-provided call-back routines, FormJacobian() and
50: FormFunction().
51: */
52: typedef struct {
53: PetscReal param; /* test problem parameter */
54: PetscInt mx; /* Discretization in x-direction */
55: PetscInt my; /* Discretization in y-direction */
56: } AppCtx;
58: /*
59: User-defined routines
60: */
61: extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
62: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
63: extern PetscErrorCode FormInitialGuess(AppCtx*,Vec);
64: extern PetscErrorCode ConvergenceTest(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*);
65: extern PetscErrorCode ConvergenceDestroy(void*);
67: int main(int argc,char **argv)
68: {
69: SNES snes; /* nonlinear solver context */
70: Vec x,r; /* solution, residual vectors */
71: Mat J; /* Jacobian matrix */
72: AppCtx user; /* user-defined application context */
74: PetscInt i,its,N,hist_its[50];
75: PetscMPIInt size;
76: PetscReal bratu_lambda_max = 6.81,bratu_lambda_min = 0.,history[50];
77: MatFDColoring fdcoloring;
78: PetscBool matrix_free = PETSC_FALSE,flg,fd_coloring = PETSC_FALSE, use_convergence_test = PETSC_FALSE;
79: KSP ksp;
80: PetscInt *testarray;
82: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
83: MPI_Comm_size(PETSC_COMM_WORLD,&size);
84: if (size != 1) SETERRQ(PETSC_COMM_SELF,1,"This is a uniprocessor example only!");
86: /*
87: Initialize problem parameters
88: */
89: user.mx = 4; user.my = 4; user.param = 6.0;
90: PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,NULL);
91: PetscOptionsGetInt(NULL,NULL,"-my",&user.my,NULL);
92: PetscOptionsGetReal(NULL,NULL,"-par",&user.param,NULL);
93: if (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min) SETERRQ(PETSC_COMM_SELF,1,"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);
105: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106: Create vector data structures; set function evaluation routine
107: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109: VecCreate(PETSC_COMM_WORLD,&x);
110: VecSetSizes(x,PETSC_DECIDE,N);
111: VecSetFromOptions(x);
112: VecDuplicate(x,&r);
114: /*
115: Set function evaluation routine and vector. Whenever the nonlinear
116: solver needs to evaluate the nonlinear function, it will call this
117: routine.
118: - Note that the final routine argument is the user-defined
119: context that provides application-specific data for the
120: function evaluation routine.
121: */
122: SNESSetFunction(snes,r,FormFunction,(void*)&user);
124: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125: Create matrix data structure; set Jacobian evaluation routine
126: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
128: /*
129: Create matrix. Here we only approximately preallocate storage space
130: for the Jacobian. See the users manual for a discussion of better
131: techniques for preallocating matrix memory.
132: */
133: PetscOptionsGetBool(NULL,NULL,"-snes_mf",&matrix_free,NULL);
134: if (!matrix_free) {
135: PetscBool matrix_free_operator = PETSC_FALSE;
136: PetscOptionsGetBool(NULL,NULL,"-snes_mf_operator",&matrix_free_operator,NULL);
137: if (matrix_free_operator) matrix_free = PETSC_FALSE;
138: }
139: if (!matrix_free) {
140: MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,5,NULL,&J);
141: }
143: /*
144: This option will cause the Jacobian to be computed via finite differences
145: efficiently using a coloring of the columns of the matrix.
146: */
147: PetscOptionsGetBool(NULL,NULL,"-snes_fd_coloring",&fd_coloring,NULL);
148: 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");
150: if (fd_coloring) {
151: ISColoring iscoloring;
152: MatColoring mc;
154: /*
155: This initializes the nonzero structure of the Jacobian. This is artificial
156: because clearly if we had a routine to compute the Jacobian we won't need
157: to use finite differences.
158: */
159: FormJacobian(snes,x,J,J,&user);
161: /*
162: Color the matrix, i.e. determine groups of columns that share no common
163: rows. These columns in the Jacobian can all be computed simulataneously.
164: */
165: MatColoringCreate(J,&mc);
166: MatColoringSetType(mc,MATCOLORINGSL);
167: MatColoringSetFromOptions(mc);
168: MatColoringApply(mc,&iscoloring);
169: MatColoringDestroy(&mc);
170: /*
171: Create the data structure that SNESComputeJacobianDefaultColor() uses
172: to compute the actual Jacobians via finite differences.
173: */
174: MatFDColoringCreate(J,iscoloring,&fdcoloring);
175: MatFDColoringSetFunction(fdcoloring,(PetscErrorCode (*)(void))FormFunction,&user);
176: MatFDColoringSetFromOptions(fdcoloring);
177: MatFDColoringSetUp(J,iscoloring,fdcoloring);
178: /*
179: Tell SNES to use the routine SNESComputeJacobianDefaultColor()
180: to compute Jacobians.
181: */
182: SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,fdcoloring);
183: ISColoringDestroy(&iscoloring);
184: }
185: /*
186: Set Jacobian matrix data structure and default Jacobian evaluation
187: routine. Whenever the nonlinear solver needs to compute the
188: Jacobian matrix, it will call this routine.
189: - Note that the final routine argument is the user-defined
190: context that provides application-specific data for the
191: Jacobian evaluation routine.
192: - The user can override with:
193: -snes_fd : default finite differencing approximation of Jacobian
194: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
195: (unless user explicitly sets preconditioner)
196: -snes_mf_operator : form preconditioning matrix as set by the user,
197: but use matrix-free approx for Jacobian-vector
198: products within Newton-Krylov method
199: */
200: else if (!matrix_free) {
201: SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);
202: }
204: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
205: Customize nonlinear solver; set runtime options
206: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
208: /*
209: Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
210: */
211: SNESSetFromOptions(snes);
213: /*
214: Set array that saves the function norms. This array is intended
215: when the user wants to save the convergence history for later use
216: rather than just to view the function norms via -snes_monitor.
217: */
218: SNESSetConvergenceHistory(snes,history,hist_its,50,PETSC_TRUE);
220: /*
221: Add a user provided convergence test; this is to test that SNESNEWTONTR properly calls the
222: user provided test before the specialized test. The convergence context is just an array to
223: test that it gets properly freed at the end
224: */
225: if (use_convergence_test) {
226: SNESGetKSP(snes,&ksp);
227: PetscMalloc1(5,&testarray);
228: KSPSetConvergenceTest(ksp,ConvergenceTest,testarray,ConvergenceDestroy);
229: }
231: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
232: Evaluate initial guess; then solve nonlinear system
233: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
234: /*
235: Note: The user should initialize the vector, x, with the initial guess
236: for the nonlinear solver prior to calling SNESSolve(). In particular,
237: to employ an initial guess of zero, the user should explicitly set
238: this vector to zero by calling VecSet().
239: */
240: FormInitialGuess(&user,x);
241: SNESSolve(snes,NULL,x);
242: SNESGetIterationNumber(snes,&its);
243: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
246: /*
247: Print the convergence history. This is intended just to demonstrate
248: use of the data attained via SNESSetConvergenceHistory().
249: */
250: PetscOptionsHasName(NULL,NULL,"-print_history",&flg);
251: if (flg) {
252: for (i=0; i<its+1; i++) {
253: PetscPrintf(PETSC_COMM_WORLD,"iteration %D: Linear iterations %D Function norm = %g\n",i,hist_its[i],(double)history[i]);
254: }
255: }
257: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
258: Free work space. All PETSc objects should be destroyed when they
259: are no longer needed.
260: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
262: if (!matrix_free) {
263: MatDestroy(&J);
264: }
265: if (fd_coloring) {
266: MatFDColoringDestroy(&fdcoloring);
267: }
268: VecDestroy(&x);
269: VecDestroy(&r);
270: SNESDestroy(&snes);
271: PetscFinalize();
272: return ierr;
273: }
274: /* ------------------------------------------------------------------- */
275: /*
276: FormInitialGuess - Forms initial approximation.
278: Input Parameters:
279: user - user-defined application context
280: X - vector
282: Output Parameter:
283: X - vector
284: */
285: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
286: {
287: 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: PetscErrorCode ierr;
343: PetscReal two = 2.0,one = 1.0,lambda,hx,hy,hxdhy,hydhx;
344: PetscScalar ut,ub,ul,ur,u,uxx,uyy,sc,*f;
345: const PetscScalar *x;
347: mx = user->mx;
348: my = user->my;
349: lambda = user->param;
350: hx = one / (PetscReal)(mx-1);
351: hy = one / (PetscReal)(my-1);
352: sc = hx*hy;
353: hxdhy = hx/hy;
354: hydhx = hy/hx;
356: /*
357: Get pointers to vector data
358: */
359: VecGetArrayRead(X,&x);
360: VecGetArray(F,&f);
362: /*
363: Compute function
364: */
365: for (j=0; j<my; j++) {
366: for (i=0; i<mx; i++) {
367: row = i + j*mx;
368: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
369: f[row] = x[row];
370: continue;
371: }
372: u = x[row];
373: ub = x[row - mx];
374: ul = x[row - 1];
375: ut = x[row + mx];
376: ur = x[row + 1];
377: uxx = (-ur + two*u - ul)*hydhx;
378: uyy = (-ut + two*u - ub)*hxdhy;
379: f[row] = uxx + uyy - sc*lambda*PetscExpScalar(u);
380: }
381: }
383: /*
384: Restore vectors
385: */
386: VecRestoreArrayRead(X,&x);
387: VecRestoreArray(F,&f);
388: return 0;
389: }
390: /* ------------------------------------------------------------------- */
391: /*
392: FormJacobian - Evaluates Jacobian matrix.
394: Input Parameters:
395: . snes - the SNES context
396: . x - input vector
397: . ptr - optional user-defined context, as set by SNESSetJacobian()
399: Output Parameters:
400: . A - Jacobian matrix
401: . B - optionally different preconditioning matrix
402: . flag - flag indicating matrix structure
403: */
404: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat jac,void *ptr)
405: {
406: AppCtx *user = (AppCtx*)ptr; /* user-defined applicatin context */
407: PetscInt i,j,row,mx,my,col[5];
408: PetscErrorCode ierr;
409: PetscScalar two = 2.0,one = 1.0,lambda,v[5],sc;
410: const PetscScalar *x;
411: PetscReal hx,hy,hxdhy,hydhx;
413: mx = user->mx;
414: my = user->my;
415: lambda = user->param;
416: hx = 1.0 / (PetscReal)(mx-1);
417: hy = 1.0 / (PetscReal)(my-1);
418: sc = hx*hy;
419: hxdhy = hx/hy;
420: hydhx = hy/hx;
422: /*
423: Get pointer to vector data
424: */
425: VecGetArrayRead(X,&x);
427: /*
428: Compute entries of the Jacobian
429: */
430: for (j=0; j<my; j++) {
431: for (i=0; i<mx; i++) {
432: row = i + j*mx;
433: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
434: MatSetValues(jac,1,&row,1,&row,&one,INSERT_VALUES);
435: continue;
436: }
437: v[0] = -hxdhy; col[0] = row - mx;
438: v[1] = -hydhx; col[1] = row - 1;
439: v[2] = two*(hydhx + hxdhy) - sc*lambda*PetscExpScalar(x[row]); col[2] = row;
440: v[3] = -hydhx; col[3] = row + 1;
441: v[4] = -hxdhy; col[4] = row + mx;
442: MatSetValues(jac,1,&row,5,col,v,INSERT_VALUES);
443: }
444: }
446: /*
447: Restore vector
448: */
449: VecRestoreArrayRead(X,&x);
451: /*
452: Assemble matrix
453: */
454: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
455: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
457: if (jac != J) {
458: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
459: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
460: }
462: return 0;
463: }
465: PetscErrorCode ConvergenceTest(KSP ksp,PetscInt it,PetscReal nrm,KSPConvergedReason *reason,void *ctx)
466: {
470: *reason = KSP_CONVERGED_ITERATING;
471: if (it > 1) {
472: *reason = KSP_CONVERGED_ITS;
473: PetscInfo(NULL,"User provided convergence test returning after 2 iterations\n");
474: }
475: return(0);
476: }
478: PetscErrorCode ConvergenceDestroy(void* ctx)
479: {
483: PetscInfo(NULL,"User provided convergence destroy called\n");
484: PetscFree(ctx);
485: return(0);
486: }
490: /*TEST
492: build:
493: requires: !single
495: test:
496: args: -ksp_gmres_cgs_refinement_type refine_always
498: test:
499: suffix: 2
500: args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always
502: test:
503: suffix: 2a
504: filter: grep -i KSPConvergedDefault > /dev/null && echo "Found KSPConvergedDefault"
505: args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always -info
506: requires: define(PETSC_USE_LOG)
508: test:
509: suffix: 2b
510: filter: grep -i "User provided convergence test" > /dev/null && echo "Found User provided convergence test"
511: args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always -use_convergence_test -info
512: requires: define(PETSC_USE_LOG)
514: test:
515: suffix: 3
516: args: -snes_monitor_short -mat_coloring_type sl -snes_fd_coloring -mx 8 -my 11 -ksp_gmres_cgs_refinement_type refine_always
518: TEST*/