Actual source code: ex5.c
1: /*$Id: ex5.c,v 1.142 2001/08/07 21:31:17 bsmith Exp $*/
3: /* Program usage: mpirun -np <procs> ex5 [-help] [all PETSc options] */
5: static char help[] = "Bratu nonlinear PDE in 2d.n
6: We solve the Bratu (SFI - solid fuel ignition) problem in a 2D rectangularn
7: domain, using distributed arrays (DAs) to partition the parallel grid.n
8: The command line options include:n
9: -par <parameter>, where <parameter> indicates the problem's nonlinearityn
10: problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81)nn";
12: /*T
13: Concepts: SNES^parallel Bratu example
14: Concepts: DA^using distributed arrays;
15: Processors: n
16: T*/
18: /* ------------------------------------------------------------------------
20: Solid Fuel Ignition (SFI) problem. This problem is modeled by
21: the partial differential equation
22:
23: -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
24:
25: with boundary conditions
26:
27: u = 0 for x = 0, x = 1, y = 0, y = 1.
28:
29: A finite difference approximation with the usual 5-point stencil
30: is used to discretize the boundary value problem to obtain a nonlinear
31: system of equations.
34: ------------------------------------------------------------------------- */
36: /*
37: Include "petscda.h" so that we can use distributed arrays (DAs).
38: Include "petscsnes.h" so that we can use SNES solvers. Note that this
39: file automatically includes:
40: petsc.h - base PETSc routines petscvec.h - vectors
41: petscsys.h - system routines petscmat.h - matrices
42: petscis.h - index sets petscksp.h - Krylov subspace methods
43: petscviewer.h - viewers petscpc.h - preconditioners
44: petscsles.h - linear solvers
45: */
46: #include petscda.h
47: #include petscsnes.h
49: /*
50: User-defined application context - contains data needed by the
51: application-provided call-back routines, FormJacobianLocal() and
52: FormFunctionLocal().
53: */
54: typedef struct {
55: DA da; /* distributed array data structure */
56: PassiveReal param; /* test problem parameter */
57: } AppCtx;
59: /*
60: User-defined routines
61: */
62: extern int FormInitialGuess(AppCtx*,Vec),FormFunctionMatlab(SNES,Vec,Vec,void*);
63: extern int FormFunctionLocal(DALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
64: extern int FormJacobianLocal(DALocalInfo*,PetscScalar**,Mat,AppCtx*);
66: #undef __FUNCT__
68: int main(int argc,char **argv)
69: {
70: SNES snes; /* nonlinear solver */
71: Vec x,r; /* solution, residual vectors */
72: Mat A,J; /* Jacobian matrix */
73: AppCtx user; /* user-defined work context */
74: int its; /* iterations for convergence */
75: PetscTruth matlab_function = PETSC_FALSE;
76: PetscTruth fd_jacobian = PETSC_FALSE,adic_jacobian=PETSC_FALSE;
77: PetscTruth adicmf_jacobian = PETSC_FALSE;
78: int ierr;
79: PetscReal bratu_lambda_max = 6.81,bratu_lambda_min = 0.;
80: MatFDColoring matfdcoloring = 0;
81: ISColoring iscoloring;
84:
85: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86: Initialize program
87: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
89: PetscInitialize(&argc,&argv,(char *)0,help);
91: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
92: Initialize problem parameters
93: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
94: user.param = 6.0;
95: PetscOptionsGetReal(PETSC_NULL,"-par",&user.param,PETSC_NULL);
96: if (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min) {
97: SETERRQ(1,"Lambda is out of range");
98: }
100: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: Create nonlinear solver context
102: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103: SNESCreate(PETSC_COMM_WORLD,&snes);
105: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106: Create distributed array (DA) to manage parallel grid and vectors
107: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
108: DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,
109: 1,1,PETSC_NULL,PETSC_NULL,&user.da);
111: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
112: Extract global vectors from DA; then duplicate for remaining
113: vectors that are the same types
114: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
115: DACreateGlobalVector(user.da,&x);
116: VecDuplicate(x,&r);
118: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119: Set local function evaluation routine
120: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
121: DASetLocalFunction(user.da,(DALocalFunction1)FormFunctionLocal);
122: DASetLocalJacobian(user.da,(DALocalFunction1)FormJacobianLocal);
123: DASetLocalAdicFunction(user.da,ad_FormFunctionLocal);
125: /* Decide which FormFunction to use */
126: PetscOptionsGetLogical(PETSC_NULL,"-matlab_function",&matlab_function,0);
128: SNESSetFunction(snes,r,SNESDAFormFunction,&user);
129: #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE)
130: if (matlab_function) {
131: SNESSetFunction(snes,r,FormFunctionMatlab,&user);
132: }
133: #endif
135: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136: Create matrix data structure; set Jacobian evaluation routine
138: Set Jacobian matrix data structure and default Jacobian evaluation
139: routine. User can override with:
140: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
141: (unless user explicitly sets preconditioner)
142: -snes_mf_operator : form preconditioning matrix as set by the user,
143: but use matrix-free approx for Jacobian-vector
144: products within Newton-Krylov method
146: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
147: DAGetMatrix(user.da,MATMPIAIJ,&J);
148: A = J;
150: PetscOptionsGetLogical(PETSC_NULL,"-fd_jacobian",&fd_jacobian,0);
151: PetscOptionsGetLogical(PETSC_NULL,"-adic_jacobian",&adic_jacobian,0);
152: PetscOptionsGetLogical(PETSC_NULL,"-adicmf_jacobian",&adicmf_jacobian,0);
153: #if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE)
154: if (adicmf_jacobian) {
155: DASetLocalAdicMFFunction(user.da,admf_FormFunctionLocal);
156: MatRegisterDAAD();
157: MatCreateDAAD(user.da,&A);
158: MatDAADSetSNES(A,snes);
159: MatDAADSetCtx(A,&user);
160: }
161: #endif
163: if (fd_jacobian) {
164: DAGetColoring(user.da,IS_COLORING_LOCAL,&iscoloring);
165: MatFDColoringCreate(J,iscoloring,&matfdcoloring);
166: ISColoringDestroy(iscoloring);
167: MatFDColoringSetFunction(matfdcoloring,(int (*)(void))SNESDAFormFunction,&user);
168: MatFDColoringSetFromOptions(matfdcoloring);
169: SNESSetJacobian(snes,A,J,SNESDefaultComputeJacobianColor,matfdcoloring);
170: #if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE)
171: } else if (adic_jacobian) {
172: DAGetColoring(user.da,IS_COLORING_GHOSTED,&iscoloring);
173: MatSetColoring(J,iscoloring);
174: ISColoringDestroy(iscoloring);
175: SNESSetJacobian(snes,A,J,SNESDAComputeJacobianWithAdic,&user);
176: #endif
177: } else {
178: SNESSetJacobian(snes,A,J,SNESDAComputeJacobian,&user);
179: }
181: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
182: Customize nonlinear solver; set runtime options
183: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
184: SNESSetFromOptions(snes);
186: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
187: Evaluate initial guess
188: Note: The user should initialize the vector, x, with the initial guess
189: for the nonlinear solver prior to calling SNESSolve(). In particular,
190: to employ an initial guess of zero, the user should explicitly set
191: this vector to zero by calling VecSet().
192: */
193: FormInitialGuess(&user,x);
195: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
196: Solve nonlinear system
197: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
198: SNESSolve(snes,x,&its);
200: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
201: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
202: PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %dn",its);
204: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
205: Free work space. All PETSc objects should be destroyed when they
206: are no longer needed.
207: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
209: if (A != J) {
210: MatDestroy(A);
211: }
212: MatDestroy(J);
213: if (matfdcoloring) {
214: MatFDColoringDestroy(matfdcoloring);
215: }
216: VecDestroy(x);
217: VecDestroy(r);
218: SNESDestroy(snes);
219: DADestroy(user.da);
220: PetscFinalize();
222: return(0);
223: }
224: /* ------------------------------------------------------------------- */
225: #undef __FUNCT__
227: /*
228: FormInitialGuess - Forms initial approximation.
230: Input Parameters:
231: user - user-defined application context
232: X - vector
234: Output Parameter:
235: X - vector
236: */
237: int FormInitialGuess(AppCtx *user,Vec X)
238: {
239: int i,j,Mx,My,ierr,xs,ys,xm,ym;
240: PetscReal lambda,temp1,temp,hx,hy;
241: PetscScalar **x;
244: DAGetInfo(user->da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
245: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
247: lambda = user->param;
248: hx = 1.0/(PetscReal)(Mx-1);
249: hy = 1.0/(PetscReal)(My-1);
250: temp1 = lambda/(lambda + 1.0);
252: /*
253: Get a pointer to vector data.
254: - For default PETSc vectors, VecGetArray() returns a pointer to
255: the data array. Otherwise, the routine is implementation dependent.
256: - You MUST call VecRestoreArray() when you no longer need access to
257: the array.
258: */
259: DAVecGetArray(user->da,X,(void**)&x);
261: /*
262: Get local grid boundaries (for 2-dimensional DA):
263: xs, ys - starting grid indices (no ghost points)
264: xm, ym - widths of local grid (no ghost points)
266: */
267: DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
269: /*
270: Compute initial guess over the locally owned part of the grid
271: */
272: for (j=ys; j<ys+ym; j++) {
273: temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
274: for (i=xs; i<xs+xm; i++) {
276: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
277: /* boundary conditions are all zero Dirichlet */
278: x[j][i] = 0.0;
279: } else {
280: x[j][i] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
281: }
282: }
283: }
285: /*
286: Restore vector
287: */
288: DAVecRestoreArray(user->da,X,(void**)&x);
290: return(0);
291: }
292: /* ------------------------------------------------------------------- */
293: #undef __FUNCT__
295: /*
296: FormFunctionLocal - Evaluates nonlinear function, F(x).
298: Process adiC: FormFunctionLocal
300: */
301: int FormFunctionLocal(DALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
302: {
303: int ierr,i,j;
304: PetscReal two = 2.0,lambda,hx,hy,hxdhy,hydhx,sc;
305: PetscScalar u,uxx,uyy;
309: lambda = user->param;
310: hx = 1.0/(PetscReal)(info->mx-1);
311: hy = 1.0/(PetscReal)(info->my-1);
312: sc = hx*hy*lambda;
313: hxdhy = hx/hy;
314: hydhx = hy/hx;
315: /*
316: Compute function over the locally owned part of the grid
317: */
318: for (j=info->ys; j<info->ys+info->ym; j++) {
319: for (i=info->xs; i<info->xs+info->xm; i++) {
320: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
321: f[j][i] = x[j][i];
322: } else {
323: u = x[j][i];
324: uxx = (two*u - x[j][i-1] - x[j][i+1])*hydhx;
325: uyy = (two*u - x[j-1][i] - x[j+1][i])*hxdhy;
326: f[j][i] = uxx + uyy - sc*PetscExpScalar(u);
327: }
328: }
329: }
331: PetscLogFlops(11*info->ym*info->xm);
332: return(0);
333: }
335: #undef __FUNCT__
337: /*
338: FormJacobianLocal - Evaluates Jacobian matrix.
341: */
342: int FormJacobianLocal(DALocalInfo *info,PetscScalar **x,Mat jac,AppCtx *user)
343: {
344: int ierr,i,j;
345: MatStencil col[5],row;
346: PetscScalar lambda,v[5],hx,hy,hxdhy,hydhx,sc;
349: lambda = user->param;
350: hx = 1.0/(PetscReal)(info->mx-1);
351: hy = 1.0/(PetscReal)(info->my-1);
352: sc = hx*hy*lambda;
353: hxdhy = hx/hy;
354: hydhx = hy/hx;
357: /*
358: Compute entries for the locally owned part of the Jacobian.
359: - Currently, all PETSc parallel matrix formats are partitioned by
360: contiguous chunks of rows across the processors.
361: - Each processor needs to insert only elements that it owns
362: locally (but any non-local elements will be sent to the
363: appropriate processor during matrix assembly).
364: - Here, we set all entries for a particular row at once.
365: - We can set matrix entries either using either
366: MatSetValuesLocal() or MatSetValues(), as discussed above.
367: */
368: for (j=info->ys; j<info->ys+info->ym; j++) {
369: for (i=info->xs; i<info->xs+info->xm; i++) {
370: row.j = j; row.i = i;
371: /* boundary points */
372: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
373: v[0] = 1.0;
374: MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
375: } else {
376: /* interior grid points */
377: v[0] = -hxdhy; col[0].j = j - 1; col[0].i = i;
378: v[1] = -hydhx; col[1].j = j; col[1].i = i-1;
379: v[2] = 2.0*(hydhx + hxdhy) - sc*PetscExpScalar(x[j][i]); col[2].j = row.j; col[2].i = row.i;
380: v[3] = -hydhx; col[3].j = j; col[3].i = i+1;
381: v[4] = -hxdhy; col[4].j = j + 1; col[4].i = i;
382: MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
383: }
384: }
385: }
387: /*
388: Assemble matrix, using the 2-step process:
389: MatAssemblyBegin(), MatAssemblyEnd().
390: */
391: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
392: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
393: /*
394: Tell the matrix we will never add a new nonzero location to the
395: matrix. If we do, it will generate an error.
396: */
397: MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR);
398: return(0);
399: }
401: /*
402: Variant of FormFunction() that computes the function in Matlab
403: */
404: #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE)
405: int FormFunctionMatlab(SNES snes,Vec X,Vec F,void *ptr)
406: {
407: AppCtx *user = (AppCtx*)ptr;
408: int ierr,Mx,My;
409: PetscReal lambda,hx,hy;
410: Vec localX,localF;
411: MPI_Comm comm;
414: DAGetLocalVector(user->da,&localX);
415: DAGetLocalVector(user->da,&localF);
416: PetscObjectSetName((PetscObject)localX,"localX");
417: PetscObjectSetName((PetscObject)localF,"localF");
418: DAGetInfo(user->da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
419: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
421: lambda = user->param;
422: hx = 1.0/(PetscReal)(Mx-1);
423: hy = 1.0/(PetscReal)(My-1);
425: PetscObjectGetComm((PetscObject)snes,&comm);
426: /*
427: Scatter ghost points to local vector,using the 2-step process
428: DAGlobalToLocalBegin(),DAGlobalToLocalEnd().
429: By placing code between these two statements, computations can be
430: done while messages are in transition.
431: */
432: DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
433: DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);
434: PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localX);
435: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(comm),"localF=ex5m(localX,%18.16e,%18.16e,%18.16e)",hx,hy,lambda);
436: PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localF);
438: /*
439: Insert values into global vector
440: */
441: DALocalToGlobal(user->da,localF,INSERT_VALUES,F);
442: DARestoreLocalVector(user->da,&localX);
443: DARestoreLocalVector(user->da,&localF);
444: return(0);
445: }
446: #endif