Actual source code: ex5.c

petsc-dev 2014-08-26
Report Typos and Errors
2: static char help[] = "Solves the modified Bratu problem in a 2D rectangular domain,\n\
3: using distributed arrays (DMDAs) to partition the parallel grid.\n\
4: The command line options include:\n\
5:   -lambda <parameter>, where <parameter> indicates the problem's nonlinearity\n\
6:   -kappa  <parameter>, where <parameter> indicates the problem's nonlinearity\n\
7:   -mx <xg>, where <xg> = number of grid points in the x-direction\n\
8:   -my <yg>, where <yg> = number of grid points in the y-direction\n\
9:   -Nx <npx>, where <npx> = number of processors in the x-direction\n\
10:   -Ny <npy>, where <npy> = number of processors in the y-direction\n\n";

12: /*T
13:    Concepts: SNES^solving a system of nonlinear equations (parallel Bratu example);
14:    Concepts: DM^using distributed arrays;
15:    Processors: n
16: T*/

18: /* ------------------------------------------------------------------------

20:     Modified Solid Fuel Ignition problem.  This problem is modeled by
21:     the partial differential equation

23:         -Laplacian u - kappa*\PartDer{u}{x} - lambda*exp(u) = 0,

25:     where

27:          0 < x,y < 1,

29:     with boundary conditions

31:              u = 0  for  x = 0, x = 1, y = 0, y = 1.

33:     A finite difference approximation with the usual 5-point stencil
34:     is used to discretize the boundary value problem to obtain a nonlinear
35:     system of equations.

37:   ------------------------------------------------------------------------- */

39: /*
40:    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
41:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
42:    file automatically includes:
43:      petscsys.h       - base PETSc routines   petscvec.h - vectors
44:      petscmat.h - matrices
45:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
46:      petscviewer.h - viewers               petscpc.h  - preconditioners
47:      petscksp.h   - linear solvers
48: */
49: #include <petscdm.h>
50: #include <petscdmda.h>
51: #include <petscsnes.h>

53: /*
54:    User-defined application context - contains data needed by the
55:    application-provided call-back routines, FormJacobian() and
56:    FormFunction().
57: */
58: typedef struct {
59:   PetscReal   param;           /* test problem parameter */
60:   PetscReal   param2;          /* test problem parameter */
61:   PetscInt    mx,my;           /* discretization in x, y directions */
62:   Vec         localX,localF;   /* ghosted local vector */
63:   DM          da;              /* distributed array data structure */
64:   PetscMPIInt rank;            /* processor rank */
65: } AppCtx;

67: /*
68:    User-defined routines
69: */
70: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*),FormInitialGuess(AppCtx*,Vec);
71: extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);

75: int main(int argc,char **argv)
76: {
77:   SNES           snes;                /* nonlinear solver */
78:   Vec            x,r;                 /* solution, residual vectors */
79:   Mat            J;                   /* Jacobian matrix */
80:   AppCtx         user;                /* user-defined work context */
81:   PetscInt       its;                 /* iterations for convergence */
82:   PetscInt       Nx,Ny;               /* number of preocessors in x- and y- directions */
83:   PetscBool      matrix_free = PETSC_FALSE;         /* flag - 1 indicates matrix-free version */
84:   PetscMPIInt    size;                /* number of processors */
85:   PetscInt       m,N;
87:   PetscReal      bratu_lambda_max = 6.81,bratu_lambda_min = 0.;
88:   PetscReal      bratu_kappa_max  = 10000,bratu_kappa_min = 0.;

90:   PetscInitialize(&argc,&argv,(char*)0,help);
91:   MPI_Comm_rank(PETSC_COMM_WORLD,&user.rank);

93:   /*
94:      Initialize problem parameters
95:   */
96:   user.mx = 4; user.my = 4; user.param = 6.0; user.param2 = 0.0;
97:   PetscOptionsGetInt(NULL,"-mx",&user.mx,NULL);
98:   PetscOptionsGetInt(NULL,"-my",&user.my,NULL);
99:   PetscOptionsGetReal(NULL,"-lambda",&user.param,NULL);
100:   if (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min) SETERRQ(PETSC_COMM_SELF,1,"Lambda is out of range");
101:   PetscOptionsGetReal(NULL,"-kappa",&user.param2,NULL);
102:   if (user.param2 >= bratu_kappa_max || user.param2 < bratu_kappa_min) SETERRQ(PETSC_COMM_SELF,1,"Kappa is out of range");
103:   PetscPrintf(PETSC_COMM_WORLD,"Solving the Bratu problem with lambda=%g, kappa=%g\n",(double)user.param,(double)user.param2);

105:   N = user.mx*user.my;

107:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108:      Create nonlinear solver context
109:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

111:   SNESCreate(PETSC_COMM_WORLD,&snes);

113:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114:      Create vector data structures; set function evaluation routine
115:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

117:   /*
118:      Create distributed array (DMDA) to manage parallel grid and vectors
119:   */
120:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
121:   Nx   = PETSC_DECIDE; Ny = PETSC_DECIDE;
122:   PetscOptionsGetInt(NULL,"-Nx",&Nx,NULL);
123:   PetscOptionsGetInt(NULL,"-Ny",&Ny,NULL);
124:   if (Nx*Ny != size && (Nx != PETSC_DECIDE || Ny != PETSC_DECIDE)) SETERRQ(PETSC_COMM_SELF,1,"Incompatible number of processors:  Nx * Ny != size");
125:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,user.mx,user.my,Nx,Ny,1,1,NULL,NULL,&user.da);
126:   SNESSetDM(snes,user.da);
127:   /*
128:      Visualize the distribution of the array across the processors
129:   */
130:   /*  DMView(user.da,PETSC_VIEWER_DRAW_WORLD); */

133:   /*
134:      Extract global and local vectors from DMDA; then duplicate for remaining
135:      vectors that are the same types
136:   */
137:   DMCreateGlobalVector(user.da,&x);
138:   DMCreateLocalVector(user.da,&user.localX);
139:   VecDuplicate(x,&r);
140:   VecDuplicate(user.localX,&user.localF);

142:   /*
143:      Set function evaluation routine and vector
144:   */
145:   SNESSetFunction(snes,r,FormFunction,(void*)&user);

147:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148:      Create matrix data structure; set Jacobian evaluation routine
149:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

151:   /*
152:      Set Jacobian matrix data structure and default Jacobian evaluation
153:      routine. User can override with:
154:      -snes_fd : default finite differencing approximation of Jacobian
155:      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
156:                 (unless user explicitly sets preconditioner)
157:      -snes_mf_operator : form preconditioning matrix as set by the user,
158:                          but use matrix-free approx for Jacobian-vector
159:                          products within Newton-Krylov method

161:      Note:  For the parallel case, vectors and matrices MUST be partitioned
162:      accordingly.  When using distributed arrays (DMDAs) to create vectors,
163:      the DMDAs determine the problem partitioning.  We must explicitly
164:      specify the local matrix dimensions upon its creation for compatibility
165:      with the vector distribution.  Thus, the generic MatCreate() routine
166:      is NOT sufficient when working with distributed arrays.

168:      Note: Here we only approximately preallocate storage space for the
169:      Jacobian.  See the users manual for a discussion of better techniques
170:      for preallocating matrix memory.
171:   */
172:   PetscOptionsGetBool(NULL,"-snes_mf",&matrix_free,NULL);
173:   if (!matrix_free) {
174:     PetscBool matrix_free_operator = PETSC_FALSE;
175:     PetscOptionsGetBool(NULL,"-snes_mf_operator",&matrix_free_operator,NULL);
176:     if (matrix_free_operator) matrix_free = PETSC_FALSE;
177:   }
178:   if (!matrix_free) {
179:     if (size == 1) {
180:       MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,5,NULL,&J);
181:     } else {
182:       VecGetLocalSize(x,&m);
183:       MatCreateAIJ(PETSC_COMM_WORLD,m,m,N,N,5,NULL,3,NULL,&J);
184:     }
185:     SNESSetJacobian(snes,J,J,FormJacobian,&user);
186:   }

188:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
189:      Customize nonlinear solver; set runtime options
190:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

192:   /*
193:      Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
194:   */
195:   SNESSetFromOptions(snes);

197:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
198:      Evaluate initial guess; then solve nonlinear system
199:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
200:   /*
201:      Note: The user should initialize the vector, x, with the initial guess
202:      for the nonlinear solver prior to calling SNESSolve().  In particular,
203:      to employ an initial guess of zero, the user should explicitly set
204:      this vector to zero by calling VecSet().
205:   */
206:   FormInitialGuess(&user,x);
207:   SNESSolve(snes,NULL,x);
208:   SNESGetIterationNumber(snes,&its);
209:   PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);

211:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
212:      Free work space.  All PETSc objects should be destroyed when they
213:      are no longer needed.
214:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

216:   if (!matrix_free) {
217:     MatDestroy(&J);
218:   }
219:   VecDestroy(&user.localX); VecDestroy(&x);
220:   VecDestroy(&user.localF); VecDestroy(&r);
221:   SNESDestroy(&snes);  DMDestroy(&user.da);
222:   PetscFinalize();

224:   return 0;
225: }
226: /* ------------------------------------------------------------------- */
229: /*
230:    FormInitialGuess - Forms initial approximation.

232:    Input Parameters:
233:    user - user-defined application context
234:    X - vector

236:    Output Parameter:
237:    X - vector
238:  */
239: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
240: {
241:   PetscInt       i,j,row,mx,my,xs,ys,xm,ym,gxm,gym,gxs,gys;
243:   PetscReal      one = 1.0,lambda,temp1,temp,hx,hy,hxdhy,hydhx,sc;
244:   PetscScalar    *x;
245:   Vec            localX = user->localX;

247:   mx    = user->mx;              my = user->my;            lambda = user->param;
248:   hx    = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1);
249:   sc    = hx*hy*lambda;       hxdhy = hx/hy;                hydhx = hy/hx;
250:   temp1 = lambda/(lambda + one);

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:   VecGetArray(localX,&x);

261:   /*
262:      Get local grid boundaries (for 2-dimensional DMDA):
263:        xs, ys   - starting grid indices (no ghost points)
264:        xm, ym   - widths of local grid (no ghost points)
265:        gxs, gys - starting grid indices (including ghost points)
266:        gxm, gym - widths of local grid (including ghost points)
267:   */
268:   DMDAGetCorners(user->da,&xs,&ys,NULL,&xm,&ym,NULL);
269:   DMDAGetGhostCorners(user->da,&gxs,&gys,NULL,&gxm,&gym,NULL);

271:   /*
272:      Compute initial guess over the locally owned part of the grid
273:   */
274:   for (j=ys; j<ys+ym; j++) {
275:     temp = (PetscReal)(PetscMin(j,my-j-1))*hy;
276:     for (i=xs; i<xs+xm; i++) {
277:       row = i - gxs + (j - gys)*gxm;
278:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
279:         x[row] = 0.0;
280:         continue;
281:       }
282:       x[row] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));
283:     }
284:   }

286:   /*
287:      Restore vector
288:   */
289:   VecRestoreArray(localX,&x);

291:   /*
292:      Insert values into global vector
293:   */
294:   DMLocalToGlobalBegin(user->da,localX,INSERT_VALUES,X);
295:   DMLocalToGlobalEnd(user->da,localX,INSERT_VALUES,X);
296:   return 0;
297: }
298: /* ------------------------------------------------------------------- */
301: /*
302:    FormFunction - Evaluates nonlinear function, F(x).

304:    Input Parameters:
305: .  snes - the SNES context
306: .  X - input vector
307: .  ptr - optional user-defined context, as set by SNESSetFunction()

309:    Output Parameter:
310: .  F - function vector
311:  */
312: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
313: {
314:   AppCtx         *user = (AppCtx*)ptr;
316:   PetscInt       i,j,row,mx,my,xs,ys,xm,ym,gxs,gys,gxm,gym;
317:   PetscReal      two = 2.0,one = 1.0,half = 0.5;
318:   PetscReal      lambda,hx,hy,hxdhy,hydhx,sc;
319:   PetscScalar    u,ux,uxx,uyy,*x,*f,kappa;
320:   Vec            localX = user->localX,localF = user->localF;

322:   mx    = user->mx;               my = user->my;        lambda = user->param;
323:   hx    = one/(PetscReal)(mx-1);  hy = one/(PetscReal)(my-1);
324:   sc    = hx*hy*lambda;        hxdhy = hx/hy;            hydhx = hy/hx;
325:   kappa = user->param2;

327:   /*
328:      Scatter ghost points to local vector, using the 2-step process
329:         DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
330:      By placing code between these two statements, computations can be
331:      done while messages are in transition.
332:   */
333:   DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
334:   DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);

336:   /*
337:      Get pointers to vector data
338:   */
339:   VecGetArray(localX,&x);
340:   VecGetArray(localF,&f);

342:   /*
343:      Get local grid boundaries
344:   */
345:   DMDAGetCorners(user->da,&xs,&ys,NULL,&xm,&ym,NULL);
346:   DMDAGetGhostCorners(user->da,&gxs,&gys,NULL,&gxm,&gym,NULL);

348:   /*
349:      Compute function over the locally owned part of the grid
350:   */
351:   for (j=ys; j<ys+ym; j++) {
352:     row = (j - gys)*gxm + xs - gxs - 1;
353:     for (i=xs; i<xs+xm; i++) {
354:       row++;
355:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
356:         f[row] = x[row];
357:         continue;
358:       }
359:       u      = x[row];
360:       ux     = (x[row+1] - x[row-1])*half*hy;
361:       uxx    = (two*u - x[row-1] - x[row+1])*hydhx;
362:       uyy    = (two*u - x[row-gxm] - x[row+gxm])*hxdhy;
363:       f[row] = uxx + uyy - kappa*ux - sc*PetscExpScalar(u);
364:     }
365:   }

367:   /*
368:      Restore vectors
369:   */
370:   VecRestoreArray(localX,&x);
371:   VecRestoreArray(localF,&f);

373:   /*
374:      Insert values into global vector
375:   */
376:   DMLocalToGlobalBegin(user->da,localF,INSERT_VALUES,F);
377:   DMLocalToGlobalEnd(user->da,localF,INSERT_VALUES,F);
378:   PetscLogFlops(11.0*ym*xm);
379:   return 0;
380: }
381: /* ------------------------------------------------------------------- */
384: /*
385:    FormJacobian - Evaluates Jacobian matrix.

387:    Input Parameters:
388: .  snes - the SNES context
389: .  x - input vector
390: .  ptr - optional user-defined context, as set by SNESSetJacobian()

392:    Output Parameters:
393: .  A - Jacobian matrix
394: .  B - optionally different preconditioning matrix
395: .  flag - flag indicating matrix structure

397:    Notes:
398:    Due to grid point reordering with DMDAs, we must always work
399:    with the local grid points, and then transform them to the new
400:    global numbering with the "ltog" mapping
401:    We cannot work directly with the global numbers for the original
402:    uniprocessor grid!
403: */
404: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat jac,void *ptr)
405: {
406:   AppCtx         *user  = (AppCtx*)ptr;   /* user-defined application context */
407:   Vec            localX = user->localX;   /* local vector */
409:   PetscInt       i,j,row,mx,my,col[5];
410:   PetscInt       xs,ys,xm,ym,gxs,gys,gxm,gym;
411:   PetscScalar    two = 2.0,one = 1.0,lambda,v[5],hx,hy,hxdhy,hydhx,sc,*x;

413:   mx = user->mx;            my = user->my;            lambda = user->param;
414:   hx = one/(PetscReal)(mx-1);  hy = one/(PetscReal)(my-1);
415:   sc = hx*hy;               hxdhy = hx/hy;            hydhx = hy/hx;

417:   /*
418:      Scatter ghost points to local vector,using the 2-step process
419:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
420:      By placing code between these two statements, computations can be
421:      done while messages are in transition.
422:   */
423:   DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
424:   DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);

426:   /*
427:      Get pointer to vector data
428:   */
429:   VecGetArray(localX,&x);

431:   /*
432:      Get local grid boundaries
433:   */
434:   DMDAGetCorners(user->da,&xs,&ys,NULL,&xm,&ym,NULL);
435:   DMDAGetGhostCorners(user->da,&gxs,&gys,NULL,&gxm,&gym,NULL);

437:   /*
438:      Compute entries for the locally owned part of the Jacobian.
439:       - Currently, all PETSc parallel matrix formats are partitioned by
440:         contiguous chunks of rows across the processors. The "grow"
441:         parameter computed below specifies the global row number
442:         corresponding to each local grid point.
443:       - Each processor needs to insert only elements that it owns
444:         locally (but any non-local elements will be sent to the
445:         appropriate processor during matrix assembly).
446:       - Always specify global row and columns of matrix entries.
447:       - Here, we set all entries for a particular row at once.
448:   */
449:   for (j=ys; j<ys+ym; j++) {
450:     row = (j - gys)*gxm + xs - gxs - 1;
451:     for (i=xs; i<xs+xm; i++) {
452:       row++;
453:       /* boundary points */
454:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
455:         MatSetValuesLocal(jac,1,&row,1,&row,&one,INSERT_VALUES);
456:         continue;
457:       }
458:       /* interior grid points */
459:       v[0] = -hxdhy; col[0] = row - gxm;
460:       v[1] = -hydhx; col[1] = row - 1;
461:       v[2] = two*(hydhx + hxdhy) - sc*lambda*PetscExpScalar(x[row]); col[2] = row;
462:       v[3] = -hydhx; col[3] = row + 1;
463:       v[4] = -hxdhy; col[4] = row + gxm;
464:       MatSetValuesLocal(jac,1,&row,5,col,v,INSERT_VALUES);
465:     }
466:   }

468:   /*
469:      Assemble matrix, using the 2-step process:
470:        MatAssemblyBegin(), MatAssemblyEnd().
471:      By placing code between these two statements, computations can be
472:      done while messages are in transition.
473:   */
474:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
475:   VecRestoreArray(localX,&x);
476:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

478:   /*
479:       Tell the matrix we will never add a new nonzero location to the
480:     matrix. If we do it will generate an error.
481:   */
482:   /* MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE); */
483:   return 0;
484: }