Actual source code: ex1.c
petsc-3.8.4 2018-03-24
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/examples/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);
63: int main(int argc,char **argv)
64: {
65: SNES snes; /* nonlinear solver context */
66: Vec x,r; /* solution, residual vectors */
67: Mat J; /* Jacobian matrix */
68: AppCtx user; /* user-defined application context */
70: PetscInt i,its,N,hist_its[50];
71: PetscMPIInt size;
72: PetscReal bratu_lambda_max = 6.81,bratu_lambda_min = 0.,history[50];
73: MatFDColoring fdcoloring;
74: PetscBool matrix_free = PETSC_FALSE,flg,fd_coloring = PETSC_FALSE;
76: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
77: MPI_Comm_size(PETSC_COMM_WORLD,&size);
78: if (size != 1) SETERRQ(PETSC_COMM_SELF,1,"This is a uniprocessor example only!");
80: /*
81: Initialize problem parameters
82: */
83: user.mx = 4; user.my = 4; user.param = 6.0;
84: PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,NULL);
85: PetscOptionsGetInt(NULL,NULL,"-my",&user.my,NULL);
86: PetscOptionsGetReal(NULL,NULL,"-par",&user.param,NULL);
87: if (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min) SETERRQ(PETSC_COMM_SELF,1,"Lambda is out of range");
88: N = user.mx*user.my;
90: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91: Create nonlinear solver context
92: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
94: SNESCreate(PETSC_COMM_WORLD,&snes);
96: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97: Create vector data structures; set function evaluation routine
98: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
100: VecCreate(PETSC_COMM_WORLD,&x);
101: VecSetSizes(x,PETSC_DECIDE,N);
102: VecSetFromOptions(x);
103: VecDuplicate(x,&r);
105: /*
106: Set function evaluation routine and vector. Whenever the nonlinear
107: solver needs to evaluate the nonlinear function, it will call this
108: routine.
109: - Note that the final routine argument is the user-defined
110: context that provides application-specific data for the
111: function evaluation routine.
112: */
113: SNESSetFunction(snes,r,FormFunction,(void*)&user);
115: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: Create matrix data structure; set Jacobian evaluation routine
117: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
119: /*
120: Create matrix. Here we only approximately preallocate storage space
121: for the Jacobian. See the users manual for a discussion of better
122: techniques for preallocating matrix memory.
123: */
124: PetscOptionsGetBool(NULL,NULL,"-snes_mf",&matrix_free,NULL);
125: if (!matrix_free) {
126: PetscBool matrix_free_operator = PETSC_FALSE;
127: PetscOptionsGetBool(NULL,NULL,"-snes_mf_operator",&matrix_free_operator,NULL);
128: if (matrix_free_operator) matrix_free = PETSC_FALSE;
129: }
130: if (!matrix_free) {
131: MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,5,NULL,&J);
132: }
134: /*
135: This option will cause the Jacobian to be computed via finite differences
136: efficiently using a coloring of the columns of the matrix.
137: */
138: PetscOptionsGetBool(NULL,NULL,"-snes_fd_coloring",&fd_coloring,NULL);
139: 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");
141: if (fd_coloring) {
142: ISColoring iscoloring;
143: MatColoring mc;
145: /*
146: This initializes the nonzero structure of the Jacobian. This is artificial
147: because clearly if we had a routine to compute the Jacobian we won't need
148: to use finite differences.
149: */
150: FormJacobian(snes,x,J,J,&user);
152: /*
153: Color the matrix, i.e. determine groups of columns that share no common
154: rows. These columns in the Jacobian can all be computed simulataneously.
155: */
156: MatColoringCreate(J,&mc);
157: MatColoringSetType(mc,MATCOLORINGSL);
158: MatColoringSetFromOptions(mc);
159: MatColoringApply(mc,&iscoloring);
160: MatColoringDestroy(&mc);
161: /*
162: Create the data structure that SNESComputeJacobianDefaultColor() uses
163: to compute the actual Jacobians via finite differences.
164: */
165: MatFDColoringCreate(J,iscoloring,&fdcoloring);
166: MatFDColoringSetFunction(fdcoloring,(PetscErrorCode (*)(void))FormFunction,&user);
167: MatFDColoringSetFromOptions(fdcoloring);
168: MatFDColoringSetUp(J,iscoloring,fdcoloring);
169: /*
170: Tell SNES to use the routine SNESComputeJacobianDefaultColor()
171: to compute Jacobians.
172: */
173: SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,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,NULL,x);
222: SNESGetIterationNumber(snes,&its);
223: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES 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(NULL,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],(double)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();
252: return ierr;
253: }
254: /* ------------------------------------------------------------------- */
255: /*
256: FormInitialGuess - Forms initial approximation.
258: Input Parameters:
259: user - user-defined application context
260: X - vector
262: Output Parameter:
263: X - vector
264: */
265: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
266: {
267: PetscInt i,j,row,mx,my;
269: PetscReal lambda,temp1,temp,hx,hy;
270: PetscScalar *x;
272: mx = user->mx;
273: my = user->my;
274: lambda = user->param;
276: hx = 1.0 / (PetscReal)(mx-1);
277: hy = 1.0 / (PetscReal)(my-1);
279: /*
280: Get a pointer to vector data.
281: - For default PETSc vectors, VecGetArray() returns a pointer to
282: the data array. Otherwise, the routine is implementation dependent.
283: - You MUST call VecRestoreArray() when you no longer need access to
284: the array.
285: */
286: VecGetArray(X,&x);
287: temp1 = lambda/(lambda + 1.0);
288: for (j=0; j<my; j++) {
289: temp = (PetscReal)(PetscMin(j,my-j-1))*hy;
290: for (i=0; i<mx; i++) {
291: row = i + j*mx;
292: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
293: x[row] = 0.0;
294: continue;
295: }
296: x[row] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));
297: }
298: }
300: /*
301: Restore vector
302: */
303: VecRestoreArray(X,&x);
304: return 0;
305: }
306: /* ------------------------------------------------------------------- */
307: /*
308: FormFunction - Evaluates nonlinear function, F(x).
310: Input Parameters:
311: . snes - the SNES context
312: . X - input vector
313: . ptr - optional user-defined context, as set by SNESSetFunction()
315: Output Parameter:
316: . F - function vector
317: */
318: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
319: {
320: AppCtx *user = (AppCtx*)ptr;
321: PetscInt i,j,row,mx,my;
322: PetscErrorCode ierr;
323: PetscReal two = 2.0,one = 1.0,lambda,hx,hy,hxdhy,hydhx;
324: PetscScalar ut,ub,ul,ur,u,uxx,uyy,sc,*f;
325: const PetscScalar *x;
327: mx = user->mx;
328: my = user->my;
329: lambda = user->param;
330: hx = one / (PetscReal)(mx-1);
331: hy = one / (PetscReal)(my-1);
332: sc = hx*hy;
333: hxdhy = hx/hy;
334: hydhx = hy/hx;
336: /*
337: Get pointers to vector data
338: */
339: VecGetArrayRead(X,&x);
340: VecGetArray(F,&f);
342: /*
343: Compute function
344: */
345: for (j=0; j<my; j++) {
346: for (i=0; i<mx; i++) {
347: row = i + j*mx;
348: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
349: f[row] = x[row];
350: continue;
351: }
352: u = x[row];
353: ub = x[row - mx];
354: ul = x[row - 1];
355: ut = x[row + mx];
356: ur = x[row + 1];
357: uxx = (-ur + two*u - ul)*hydhx;
358: uyy = (-ut + two*u - ub)*hxdhy;
359: f[row] = uxx + uyy - sc*lambda*PetscExpScalar(u);
360: }
361: }
363: /*
364: Restore vectors
365: */
366: VecRestoreArrayRead(X,&x);
367: VecRestoreArray(F,&f);
368: return 0;
369: }
370: /* ------------------------------------------------------------------- */
371: /*
372: FormJacobian - Evaluates Jacobian matrix.
374: Input Parameters:
375: . snes - the SNES context
376: . x - input vector
377: . ptr - optional user-defined context, as set by SNESSetJacobian()
379: Output Parameters:
380: . A - Jacobian matrix
381: . B - optionally different preconditioning matrix
382: . flag - flag indicating matrix structure
383: */
384: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat jac,void *ptr)
385: {
386: AppCtx *user = (AppCtx*)ptr; /* user-defined applicatin context */
387: PetscInt i,j,row,mx,my,col[5];
388: PetscErrorCode ierr;
389: PetscScalar two = 2.0,one = 1.0,lambda,v[5],sc;
390: const PetscScalar *x;
391: PetscReal hx,hy,hxdhy,hydhx;
393: mx = user->mx;
394: my = user->my;
395: lambda = user->param;
396: hx = 1.0 / (PetscReal)(mx-1);
397: hy = 1.0 / (PetscReal)(my-1);
398: sc = hx*hy;
399: hxdhy = hx/hy;
400: hydhx = hy/hx;
402: /*
403: Get pointer to vector data
404: */
405: VecGetArrayRead(X,&x);
407: /*
408: Compute entries of the Jacobian
409: */
410: for (j=0; j<my; j++) {
411: for (i=0; i<mx; i++) {
412: row = i + j*mx;
413: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
414: MatSetValues(jac,1,&row,1,&row,&one,INSERT_VALUES);
415: continue;
416: }
417: v[0] = -hxdhy; col[0] = row - mx;
418: v[1] = -hydhx; col[1] = row - 1;
419: v[2] = two*(hydhx + hxdhy) - sc*lambda*PetscExpScalar(x[row]); col[2] = row;
420: v[3] = -hydhx; col[3] = row + 1;
421: v[4] = -hxdhy; col[4] = row + mx;
422: MatSetValues(jac,1,&row,5,col,v,INSERT_VALUES);
423: }
424: }
426: /*
427: Restore vector
428: */
429: VecRestoreArrayRead(X,&x);
431: /*
432: Assemble matrix
433: */
434: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
435: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
437: if (jac != J) {
438: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
439: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
440: }
442: return 0;
443: }