Actual source code: ex5.c
2: static char help[] = "Bratu nonlinear PDE in 2d.\n\
3: We solve the Bratu (SFI - solid fuel ignition) problem in a 2D rectangular\n\
4: domain, using distributed arrays (DAs) to partition the parallel grid.\n\
5: The command line options include:\n\
6: -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
7: problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81)\n\n";
9: /*T
10: Concepts: SNES^parallel Bratu example
11: Concepts: DA^using distributed arrays;
12: Concepts: IS coloirng types;
13: Processors: n
14: T*/
16: /* ------------------------------------------------------------------------
18: Solid Fuel Ignition (SFI) problem. This problem is modeled by
19: the partial differential equation
20:
21: -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
22:
23: with boundary conditions
24:
25: u = 0 for x = 0, x = 1, y = 0, y = 1.
26:
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: Program usage: mpiexec -n <procs> ex5 [-help] [all PETSc options]
32: e.g.,
33: ./ex5 -fd_jacobian -mat_fd_coloring_view_draw -draw_pause -1
34: mpiexec -n 2 ./ex5 -fd_jacobian_ghosted -log_summary
36: ------------------------------------------------------------------------- */
38: /*
39: Include "petscda.h" so that we can use distributed arrays (DAs).
40: Include "petscsnes.h" so that we can use SNES solvers. Note that this
41: file automatically includes:
42: petscsys.h - base PETSc routines petscvec.h - vectors
43: petscmat.h - matrices
44: petscis.h - index sets petscksp.h - Krylov subspace methods
45: petscviewer.h - viewers petscpc.h - preconditioners
46: petscksp.h - linear solvers
47: */
48: #include petscda.h
49: #include petscsnes.h
51: /*
52: User-defined application context - contains data needed by the
53: application-provided call-back routines, FormJacobianLocal() and
54: FormFunctionLocal().
55: */
56: typedef struct {
57: DA da; /* distributed array data structure */
58: PassiveReal param; /* test problem parameter */
59: } AppCtx;
61: /*
62: User-defined routines
63: */
71: int main(int argc,char **argv)
72: {
73: SNES snes; /* nonlinear solver */
74: Vec x,r; /* solution, residual vectors */
75: Mat A,J; /* Jacobian matrix */
76: AppCtx user; /* user-defined work context */
77: PetscInt its; /* iterations for convergence */
78: PetscTruth matlab_function = PETSC_FALSE;
79: PetscTruth fd_jacobian = PETSC_FALSE,adic_jacobian=PETSC_FALSE,fd_jacobian_ghosted=PETSC_FALSE;
80: PetscTruth adicmf_jacobian = PETSC_FALSE;
81: PetscErrorCode ierr;
82: PetscReal bratu_lambda_max = 6.81,bratu_lambda_min = 0.;
83: MatFDColoring matfdcoloring = 0;
84: ISColoring iscoloring;
86: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
87: Initialize program
88: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
90: PetscInitialize(&argc,&argv,(char *)0,help);
92: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
93: Initialize problem parameters
94: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
95: user.param = 6.0;
96: PetscOptionsGetReal(PETSC_NULL,"-par",&user.param,PETSC_NULL);
97: if (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min) {
98: SETERRQ(1,"Lambda is out of range");
99: }
101: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102: Create nonlinear solver context
103: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
104: SNESCreate(PETSC_COMM_WORLD,&snes);
106: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: Create distributed array (DA) to manage parallel grid and vectors
108: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109: DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,
110: 1,1,PETSC_NULL,PETSC_NULL,&user.da);
111: DASetUniformCoordinates(user.da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
113: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: Extract global vectors from DA; then duplicate for remaining
115: vectors that are the same types
116: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
117: DACreateGlobalVector(user.da,&x);
118: VecDuplicate(x,&r);
120: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121: Create matrix data structure; set Jacobian evaluation routine
123: Set Jacobian matrix data structure and default Jacobian evaluation
124: routine. User can override with:
125: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
126: (unless user explicitly sets preconditioner)
127: -snes_mf_operator : form preconditioning matrix as set by the user,
128: but use matrix-free approx for Jacobian-vector
129: products within Newton-Krylov method
131: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
132: /* J can be type of MATAIJ,MATBAIJ or MATSBAIJ */
133: DAGetMatrix(user.da,MATAIJ,&J);
134:
135: A = J;
136: PetscOptionsGetTruth(PETSC_NULL,"-fd_jacobian",&fd_jacobian,0);
137: PetscOptionsGetTruth(PETSC_NULL,"-fd_jacobian_ghosted",&fd_jacobian_ghosted,0);
138: PetscOptionsGetTruth(PETSC_NULL,"-adic_jacobian",&adic_jacobian,0);
139: PetscOptionsGetTruth(PETSC_NULL,"-adicmf_jacobian",&adicmf_jacobian,0);
140: #if defined(PETSC_HAVE_ADIC)
141: if (adicmf_jacobian) {
142: DASetLocalAdicMFFunction(user.da,admf_FormFunctionLocal);
143: MatRegisterDAAD();
144: MatCreateDAAD(user.da,&A);
145: MatDAADSetSNES(A,snes);
146: MatDAADSetCtx(A,&user);
147: }
148: #endif
150: if (fd_jacobian) {
151: DAGetColoring(user.da,IS_COLORING_GLOBAL,MATAIJ,&iscoloring);
152: MatFDColoringCreate(J,iscoloring,&matfdcoloring);
153: ISColoringDestroy(iscoloring);
154: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))SNESDAFormFunction,&user);
155: MatFDColoringSetFromOptions(matfdcoloring);
156: SNESSetJacobian(snes,A,J,SNESDefaultComputeJacobianColor,matfdcoloring);
157: } else if (fd_jacobian_ghosted) {
158: DAGetColoring(user.da,IS_COLORING_GHOSTED,MATAIJ,&iscoloring);
159: MatFDColoringCreate(J,iscoloring,&matfdcoloring);
160: ISColoringDestroy(iscoloring);
161: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))SNESDAFormFunction,&user);
162: MatFDColoringSetFromOptions(matfdcoloring);
163: /* now, call a customized SNESDefaultComputeJacobianColor() */
164: SNESSetJacobian(snes,A,J,MySNESDefaultComputeJacobianColor,matfdcoloring);
165: #if defined(PETSC_HAVE_ADIC)
166: } else if (adic_jacobian) {
167: DAGetColoring(user.da,IS_COLORING_GHOSTED,MATAIJ,&iscoloring);
168: MatSetColoring(J,iscoloring);
169: ISColoringDestroy(iscoloring);
170: SNESSetJacobian(snes,A,J,SNESDAComputeJacobianWithAdic,&user);
171: #endif
172: } else {
173: SNESSetJacobian(snes,A,J,SNESDAComputeJacobian,&user);
174: }
176: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
177: Set local function evaluation routine
178: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
179: DASetLocalFunction(user.da,(DALocalFunction1)FormFunctionLocal);
180: DASetLocalJacobian(user.da,(DALocalFunction1)FormJacobianLocal);
181: DASetLocalAdicFunction(user.da,ad_FormFunctionLocal);
183: /* Decide which FormFunction to use */
184: PetscOptionsGetTruth(PETSC_NULL,"-matlab_function",&matlab_function,0);
186: SNESSetFunction(snes,r,SNESDAFormFunction,&user);
187: #if defined(PETSC_HAVE_MATLAB_ENGINE)
188: if (matlab_function) {
189: SNESSetFunction(snes,r,FormFunctionMatlab,&user);
190: }
191: #endif
193: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
194: Customize nonlinear solver; set runtime options
195: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
196: SNESSetFromOptions(snes);
198: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
199: Evaluate initial guess
200: Note: The user should initialize the vector, x, with the initial guess
201: for the nonlinear solver prior to calling SNESSolve(). In particular,
202: to employ an initial guess of zero, the user should explicitly set
203: this vector to zero by calling VecSet().
204: */
206: {
207: PetscTruth test_appctx = PETSC_FALSE;
208: PetscOptionsGetTruth(PETSC_NULL,"-test_appctx",&test_appctx,0);
209: if (test_appctx) {
210: AppCtx *puser;
211: SNESSetApplicationContext(snes,&user);
212: SNESGetApplicationContext(snes,(void **)&puser);
213: FormInitialGuess(puser,x);
214: } else {
215: FormInitialGuess(&user,x);
216: }
217: }
219: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
220: Solve nonlinear system
221: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
222: SNESSolve(snes,PETSC_NULL,x);
223: SNESGetIterationNumber(snes,&its);
225: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
226: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
227: PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %D\n",its);
228: {
229: PetscViewer viewer;
231: PetscViewerCreate(PETSC_COMM_WORLD, &viewer);
232: //PetscViewerSetType(viewer, PETSC_VIEWER_DRAW);
233: //PetscViewerDrawSetInfo(viewer, PETSC_NULL, "Solution", PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE);
234: PetscViewerSetType(viewer, PETSC_VIEWER_ASCII);
235: PetscViewerFileSetName(viewer, "ex5_sol.vtk");
236: PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_VTK);
237: DAView(user.da, viewer);
238: VecView(x, viewer);
239: PetscViewerDestroy(viewer);
240: }
242: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
243: Free work space. All PETSc objects should be destroyed when they
244: are no longer needed.
245: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
247: if (A != J) {
248: MatDestroy(A);
249: }
250: MatDestroy(J);
251: if (matfdcoloring) {
252: MatFDColoringDestroy(matfdcoloring);
253: }
254: VecDestroy(x);
255: VecDestroy(r);
256: SNESDestroy(snes);
257: DADestroy(user.da);
258: PetscFinalize();
260: return(0);
261: }
262: /* ------------------------------------------------------------------- */
265: /*
266: FormInitialGuess - Forms initial approximation.
268: Input Parameters:
269: user - user-defined application context
270: X - vector
272: Output Parameter:
273: X - vector
274: */
275: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
276: {
277: PetscInt i,j,Mx,My,xs,ys,xm,ym;
279: PetscReal lambda,temp1,temp,hx,hy;
280: PetscScalar **x;
283: DAGetInfo(user->da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
284: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
286: lambda = user->param;
287: hx = 1.0/(PetscReal)(Mx-1);
288: hy = 1.0/(PetscReal)(My-1);
289: temp1 = lambda/(lambda + 1.0);
291: /*
292: Get a pointer to vector data.
293: - For default PETSc vectors, VecGetArray() returns a pointer to
294: the data array. Otherwise, the routine is implementation dependent.
295: - You MUST call VecRestoreArray() when you no longer need access to
296: the array.
297: */
298: DAVecGetArray(user->da,X,&x);
300: /*
301: Get local grid boundaries (for 2-dimensional DA):
302: xs, ys - starting grid indices (no ghost points)
303: xm, ym - widths of local grid (no ghost points)
305: */
306: DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
308: /*
309: Compute initial guess over the locally owned part of the grid
310: */
311: for (j=ys; j<ys+ym; j++) {
312: temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
313: for (i=xs; i<xs+xm; i++) {
314: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
315: /* boundary conditions are all zero Dirichlet */
316: x[j][i] = 0.0;
317: } else {
318: x[j][i] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
319: }
320: }
321: }
323: /*
324: Restore vector
325: */
326: DAVecRestoreArray(user->da,X,&x);
328: return(0);
329: }
330: /* ------------------------------------------------------------------- */
333: /*
334: FormFunctionLocal - Evaluates nonlinear function, F(x).
336: Process adiC(36): FormFunctionLocal
338: */
339: PetscErrorCode FormFunctionLocal(DALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
340: {
342: PetscInt i,j;
343: PetscReal lambda,hx,hy,hxdhy,hydhx,sc;
344: PetscScalar u,uxx,uyy;
348: lambda = user->param;
349: hx = 1.0/(PetscReal)(info->mx-1);
350: hy = 1.0/(PetscReal)(info->my-1);
351: sc = hx*hy*lambda;
352: hxdhy = hx/hy;
353: hydhx = hy/hx;
354: /*
355: Compute function over the locally owned part of the grid
356: */
357: for (j=info->ys; j<info->ys+info->ym; j++) {
358: for (i=info->xs; i<info->xs+info->xm; i++) {
359: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
360: f[j][i] = x[j][i];
361: } else {
362: u = x[j][i];
363: uxx = (2.0*u - x[j][i-1] - x[j][i+1])*hydhx;
364: uyy = (2.0*u - x[j-1][i] - x[j+1][i])*hxdhy;
365: f[j][i] = uxx + uyy - sc*PetscExpScalar(u);
366: }
367: }
368: }
370: PetscLogFlops(11.0*info->ym*info->xm);
371: return(0);
372: }
376: /*
377: FormJacobianLocal - Evaluates Jacobian matrix.
378: */
379: PetscErrorCode FormJacobianLocal(DALocalInfo *info,PetscScalar **x,Mat jac,AppCtx *user)
380: {
382: PetscInt i,j;
383: MatStencil col[5],row;
384: PetscScalar lambda,v[5],hx,hy,hxdhy,hydhx,sc;
387: lambda = user->param;
388: hx = 1.0/(PetscReal)(info->mx-1);
389: hy = 1.0/(PetscReal)(info->my-1);
390: sc = hx*hy*lambda;
391: hxdhy = hx/hy;
392: hydhx = hy/hx;
395: /*
396: Compute entries for the locally owned part of the Jacobian.
397: - Currently, all PETSc parallel matrix formats are partitioned by
398: contiguous chunks of rows across the processors.
399: - Each processor needs to insert only elements that it owns
400: locally (but any non-local elements will be sent to the
401: appropriate processor during matrix assembly).
402: - Here, we set all entries for a particular row at once.
403: - We can set matrix entries either using either
404: MatSetValuesLocal() or MatSetValues(), as discussed above.
405: */
406: for (j=info->ys; j<info->ys+info->ym; j++) {
407: for (i=info->xs; i<info->xs+info->xm; i++) {
408: row.j = j; row.i = i;
409: /* boundary points */
410: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
411: v[0] = 1.0;
412: MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
413: } else {
414: /* interior grid points */
415: v[0] = -hxdhy; col[0].j = j - 1; col[0].i = i;
416: v[1] = -hydhx; col[1].j = j; col[1].i = i-1;
417: v[2] = 2.0*(hydhx + hxdhy) - sc*PetscExpScalar(x[j][i]); col[2].j = row.j; col[2].i = row.i;
418: v[3] = -hydhx; col[3].j = j; col[3].i = i+1;
419: v[4] = -hxdhy; col[4].j = j + 1; col[4].i = i;
420: MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
421: }
422: }
423: }
425: /*
426: Assemble matrix, using the 2-step process:
427: MatAssemblyBegin(), MatAssemblyEnd().
428: */
429: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
430: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
431: /*
432: Tell the matrix we will never add a new nonzero location to the
433: matrix. If we do, it will generate an error.
434: */
435: MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
436: return(0);
437: }
439: /*
440: Variant of FormFunction() that computes the function in Matlab
441: */
442: #if defined(PETSC_HAVE_MATLAB_ENGINE)
443: PetscErrorCode FormFunctionMatlab(SNES snes,Vec X,Vec F,void *ptr)
444: {
445: AppCtx *user = (AppCtx*)ptr;
447: PetscInt Mx,My;
448: PetscReal lambda,hx,hy;
449: Vec localX,localF;
450: MPI_Comm comm;
453: DAGetLocalVector(user->da,&localX);
454: DAGetLocalVector(user->da,&localF);
455: PetscObjectSetName((PetscObject)localX,"localX");
456: PetscObjectSetName((PetscObject)localF,"localF");
457: DAGetInfo(user->da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
458: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
460: lambda = user->param;
461: hx = 1.0/(PetscReal)(Mx-1);
462: hy = 1.0/(PetscReal)(My-1);
464: PetscObjectGetComm((PetscObject)snes,&comm);
465: /*
466: Scatter ghost points to local vector,using the 2-step process
467: DAGlobalToLocalBegin(),DAGlobalToLocalEnd().
468: By placing code between these two statements, computations can be
469: done while messages are in transition.
470: */
471: DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
472: DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);
473: PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localX);
474: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(comm),"localF=ex5m(localX,%18.16e,%18.16e,%18.16e)",hx,hy,lambda);
475: PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localF);
477: /*
478: Insert values into global vector
479: */
480: DALocalToGlobal(user->da,localF,INSERT_VALUES,F);
481: DARestoreLocalVector(user->da,&localX);
482: DARestoreLocalVector(user->da,&localF);
483: return(0);
484: }
485: #endif
489: /*
490: MySNESDefaultComputeJacobianColor - Computes the Jacobian using
491: finite differences and coloring to exploit matrix sparsity.
492: It is customized from SNESDefaultComputeJacobianColor.
493: The input global vector x1 is scattered to x1_local
494: which then is passed into MatFDColoringApply() for reducing the
495: VecScatterBingin/End.
496: */
497: PetscErrorCode MySNESDefaultComputeJacobianColor(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx)
498: {
499: MatFDColoring color = (MatFDColoring) ctx;
501: Vec f;
502: PetscErrorCode (*ff)(void),(*fd)(void);
503: void *fctx;
504: DA da;
505: Vec x1_loc;
508: *flag = SAME_NONZERO_PATTERN;
509: SNESGetFunction(snes,&f,(PetscErrorCode (**)(SNES,Vec,Vec,void*))&ff,0);
510: MatFDColoringGetFunction(color,&fd,&fctx);
511: if (fd == ff) { /* reuse function value computed in SNES */
512: MatFDColoringSetF(color,f);
513: }
514: /* Now, get x1_loc and scatter global x1 onto x1_loc */
515: da = *(DA*)fctx;
516: DAGetLocalVector(da,&x1_loc);
517: DAGlobalToLocalBegin(da,x1,INSERT_VALUES,x1_loc);
518: DAGlobalToLocalEnd(da,x1,INSERT_VALUES,x1_loc);
519: MatFDColoringApply(*B,color,x1_loc,flag,snes);
520: DARestoreLocalVector(da,&x1_loc);
521: if (*J != *B) {
522: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
523: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
524: }
525: return(0);
526: }