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