Actual source code: ex1.c
petsc-3.9.4 2018-09-11
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);
65: int main(int argc,char **argv)
66: {
67: SNES snes; /* nonlinear solver context */
68: Vec x,r; /* solution, residual vectors */
69: Mat J; /* Jacobian matrix */
70: 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;
78: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
79: MPI_Comm_size(PETSC_COMM_WORLD,&size);
80: if (size != 1) SETERRQ(PETSC_COMM_SELF,1,"This is a uniprocessor example only!");
82: /*
83: Initialize problem parameters
84: */
85: user.mx = 4; user.my = 4; user.param = 6.0;
86: PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,NULL);
87: PetscOptionsGetInt(NULL,NULL,"-my",&user.my,NULL);
88: PetscOptionsGetReal(NULL,NULL,"-par",&user.param,NULL);
89: if (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min) SETERRQ(PETSC_COMM_SELF,1,"Lambda is out of range");
90: N = user.mx*user.my;
92: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
93: Create nonlinear solver context
94: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
96: SNESCreate(PETSC_COMM_WORLD,&snes);
98: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
99: Create vector data structures; set function evaluation routine
100: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
102: VecCreate(PETSC_COMM_WORLD,&x);
103: VecSetSizes(x,PETSC_DECIDE,N);
104: VecSetFromOptions(x);
105: VecDuplicate(x,&r);
107: /*
108: Set function evaluation routine and vector. Whenever the nonlinear
109: solver needs to evaluate the nonlinear function, it will call this
110: routine.
111: - Note that the final routine argument is the user-defined
112: context that provides application-specific data for the
113: function evaluation routine.
114: */
115: SNESSetFunction(snes,r,FormFunction,(void*)&user);
117: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118: Create matrix data structure; set Jacobian evaluation routine
119: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
121: /*
122: Create matrix. Here we only approximately preallocate storage space
123: for the Jacobian. See the users manual for a discussion of better
124: techniques for preallocating matrix memory.
125: */
126: PetscOptionsGetBool(NULL,NULL,"-snes_mf",&matrix_free,NULL);
127: if (!matrix_free) {
128: PetscBool matrix_free_operator = PETSC_FALSE;
129: PetscOptionsGetBool(NULL,NULL,"-snes_mf_operator",&matrix_free_operator,NULL);
130: if (matrix_free_operator) matrix_free = PETSC_FALSE;
131: }
132: if (!matrix_free) {
133: MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,5,NULL,&J);
134: }
136: /*
137: This option will cause the Jacobian to be computed via finite differences
138: efficiently using a coloring of the columns of the matrix.
139: */
140: PetscOptionsGetBool(NULL,NULL,"-snes_fd_coloring",&fd_coloring,NULL);
141: 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");
143: if (fd_coloring) {
144: ISColoring iscoloring;
145: MatColoring mc;
147: /*
148: This initializes the nonzero structure of the Jacobian. This is artificial
149: because clearly if we had a routine to compute the Jacobian we won't need
150: to use finite differences.
151: */
152: FormJacobian(snes,x,J,J,&user);
154: /*
155: Color the matrix, i.e. determine groups of columns that share no common
156: rows. These columns in the Jacobian can all be computed simulataneously.
157: */
158: MatColoringCreate(J,&mc);
159: MatColoringSetType(mc,MATCOLORINGSL);
160: MatColoringSetFromOptions(mc);
161: MatColoringApply(mc,&iscoloring);
162: MatColoringDestroy(&mc);
163: /*
164: Create the data structure that SNESComputeJacobianDefaultColor() uses
165: to compute the actual Jacobians via finite differences.
166: */
167: MatFDColoringCreate(J,iscoloring,&fdcoloring);
168: MatFDColoringSetFunction(fdcoloring,(PetscErrorCode (*)(void))FormFunction,&user);
169: MatFDColoringSetFromOptions(fdcoloring);
170: MatFDColoringSetUp(J,iscoloring,fdcoloring);
171: /*
172: Tell SNES to use the routine SNESComputeJacobianDefaultColor()
173: to compute Jacobians.
174: */
175: SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,fdcoloring);
176: ISColoringDestroy(&iscoloring);
177: }
178: /*
179: Set Jacobian matrix data structure and default Jacobian evaluation
180: routine. Whenever the nonlinear solver needs to compute the
181: Jacobian matrix, it will call this routine.
182: - Note that the final routine argument is the user-defined
183: context that provides application-specific data for the
184: Jacobian evaluation routine.
185: - The user can override with:
186: -snes_fd : default finite differencing approximation of Jacobian
187: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
188: (unless user explicitly sets preconditioner)
189: -snes_mf_operator : form preconditioning matrix as set by the user,
190: but use matrix-free approx for Jacobian-vector
191: products within Newton-Krylov method
192: */
193: else if (!matrix_free) {
194: SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);
195: }
197: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
198: Customize nonlinear solver; set runtime options
199: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
201: /*
202: Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
203: */
204: SNESSetFromOptions(snes);
206: /*
207: Set array that saves the function norms. This array is intended
208: when the user wants to save the convergence history for later use
209: rather than just to view the function norms via -snes_monitor.
210: */
211: SNESSetConvergenceHistory(snes,history,hist_its,50,PETSC_TRUE);
213: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
214: Evaluate initial guess; then solve nonlinear system
215: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
216: /*
217: Note: The user should initialize the vector, x, with the initial guess
218: for the nonlinear solver prior to calling SNESSolve(). In particular,
219: to employ an initial guess of zero, the user should explicitly set
220: this vector to zero by calling VecSet().
221: */
222: FormInitialGuess(&user,x);
223: SNESSolve(snes,NULL,x);
224: SNESGetIterationNumber(snes,&its);
225: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
228: /*
229: Print the convergence history. This is intended just to demonstrate
230: use of the data attained via SNESSetConvergenceHistory().
231: */
232: PetscOptionsHasName(NULL,NULL,"-print_history",&flg);
233: if (flg) {
234: for (i=0; i<its+1; i++) {
235: PetscPrintf(PETSC_COMM_WORLD,"iteration %D: Linear iterations %D Function norm = %g\n",i,hist_its[i],(double)history[i]);
236: }
237: }
239: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
240: Free work space. All PETSc objects should be destroyed when they
241: are no longer needed.
242: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
244: if (!matrix_free) {
245: MatDestroy(&J);
246: }
247: if (fd_coloring) {
248: MatFDColoringDestroy(&fdcoloring);
249: }
250: VecDestroy(&x);
251: VecDestroy(&r);
252: SNESDestroy(&snes);
253: PetscFinalize();
254: return ierr;
255: }
256: /* ------------------------------------------------------------------- */
257: /*
258: FormInitialGuess - Forms initial approximation.
260: Input Parameters:
261: user - user-defined application context
262: X - vector
264: Output Parameter:
265: X - vector
266: */
267: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
268: {
269: PetscInt i,j,row,mx,my;
271: PetscReal lambda,temp1,temp,hx,hy;
272: PetscScalar *x;
274: mx = user->mx;
275: my = user->my;
276: lambda = user->param;
278: hx = 1.0 / (PetscReal)(mx-1);
279: hy = 1.0 / (PetscReal)(my-1);
281: /*
282: Get a pointer to vector data.
283: - For default PETSc vectors, VecGetArray() returns a pointer to
284: the data array. Otherwise, the routine is implementation dependent.
285: - You MUST call VecRestoreArray() when you no longer need access to
286: the array.
287: */
288: VecGetArray(X,&x);
289: temp1 = lambda/(lambda + 1.0);
290: for (j=0; j<my; j++) {
291: temp = (PetscReal)(PetscMin(j,my-j-1))*hy;
292: for (i=0; i<mx; i++) {
293: row = i + j*mx;
294: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
295: x[row] = 0.0;
296: continue;
297: }
298: x[row] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));
299: }
300: }
302: /*
303: Restore vector
304: */
305: VecRestoreArray(X,&x);
306: return 0;
307: }
308: /* ------------------------------------------------------------------- */
309: /*
310: FormFunction - Evaluates nonlinear function, F(x).
312: Input Parameters:
313: . snes - the SNES context
314: . X - input vector
315: . ptr - optional user-defined context, as set by SNESSetFunction()
317: Output Parameter:
318: . F - function vector
319: */
320: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
321: {
322: AppCtx *user = (AppCtx*)ptr;
323: PetscInt i,j,row,mx,my;
324: PetscErrorCode ierr;
325: PetscReal two = 2.0,one = 1.0,lambda,hx,hy,hxdhy,hydhx;
326: PetscScalar ut,ub,ul,ur,u,uxx,uyy,sc,*f;
327: const PetscScalar *x;
329: mx = user->mx;
330: my = user->my;
331: lambda = user->param;
332: hx = one / (PetscReal)(mx-1);
333: hy = one / (PetscReal)(my-1);
334: sc = hx*hy;
335: hxdhy = hx/hy;
336: hydhx = hy/hx;
338: /*
339: Get pointers to vector data
340: */
341: VecGetArrayRead(X,&x);
342: VecGetArray(F,&f);
344: /*
345: Compute function
346: */
347: for (j=0; j<my; j++) {
348: for (i=0; i<mx; i++) {
349: row = i + j*mx;
350: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
351: f[row] = x[row];
352: continue;
353: }
354: u = x[row];
355: ub = x[row - mx];
356: ul = x[row - 1];
357: ut = x[row + mx];
358: ur = x[row + 1];
359: uxx = (-ur + two*u - ul)*hydhx;
360: uyy = (-ut + two*u - ub)*hxdhy;
361: f[row] = uxx + uyy - sc*lambda*PetscExpScalar(u);
362: }
363: }
365: /*
366: Restore vectors
367: */
368: VecRestoreArrayRead(X,&x);
369: VecRestoreArray(F,&f);
370: return 0;
371: }
372: /* ------------------------------------------------------------------- */
373: /*
374: FormJacobian - Evaluates Jacobian matrix.
376: Input Parameters:
377: . snes - the SNES context
378: . x - input vector
379: . ptr - optional user-defined context, as set by SNESSetJacobian()
381: Output Parameters:
382: . A - Jacobian matrix
383: . B - optionally different preconditioning matrix
384: . flag - flag indicating matrix structure
385: */
386: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat jac,void *ptr)
387: {
388: AppCtx *user = (AppCtx*)ptr; /* user-defined applicatin context */
389: PetscInt i,j,row,mx,my,col[5];
390: PetscErrorCode ierr;
391: PetscScalar two = 2.0,one = 1.0,lambda,v[5],sc;
392: const PetscScalar *x;
393: PetscReal hx,hy,hxdhy,hydhx;
395: mx = user->mx;
396: my = user->my;
397: lambda = user->param;
398: hx = 1.0 / (PetscReal)(mx-1);
399: hy = 1.0 / (PetscReal)(my-1);
400: sc = hx*hy;
401: hxdhy = hx/hy;
402: hydhx = hy/hx;
404: /*
405: Get pointer to vector data
406: */
407: VecGetArrayRead(X,&x);
409: /*
410: Compute entries of the Jacobian
411: */
412: for (j=0; j<my; j++) {
413: for (i=0; i<mx; i++) {
414: row = i + j*mx;
415: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
416: MatSetValues(jac,1,&row,1,&row,&one,INSERT_VALUES);
417: continue;
418: }
419: v[0] = -hxdhy; col[0] = row - mx;
420: v[1] = -hydhx; col[1] = row - 1;
421: v[2] = two*(hydhx + hxdhy) - sc*lambda*PetscExpScalar(x[row]); col[2] = row;
422: v[3] = -hydhx; col[3] = row + 1;
423: v[4] = -hxdhy; col[4] = row + mx;
424: MatSetValues(jac,1,&row,5,col,v,INSERT_VALUES);
425: }
426: }
428: /*
429: Restore vector
430: */
431: VecRestoreArrayRead(X,&x);
433: /*
434: Assemble matrix
435: */
436: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
437: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
439: if (jac != J) {
440: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
441: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
442: }
444: return 0;
445: }
449: /*TEST
451: build:
452: requires: !single
454: test:
455: args: -ksp_gmres_cgs_refinement_type refine_always
457: test:
458: suffix: 2
459: args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always
461: test:
462: suffix: 3
463: args: -snes_monitor_short -mat_coloring_type sl -snes_fd_coloring -mx 8 -my 11 -ksp_gmres_cgs_refinement_type refine_always
465: TEST*/