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 rectangular\n\
  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 nonlinearity\n\
 10:      problem SFI:  <parameter> = Bratu parameter (0 <= par <= 6.81)\n\n";

 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:      petscksp.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*);

 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) && !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,MATAIJ,&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);
199:   SNESGetIterationNumber(snes,&its);

201:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
203:   PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %d\n",its);

205:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
206:      Free work space.  All PETSc objects should be destroyed when they
207:      are no longer needed.
208:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

210:   if (A != J) {
211:     MatDestroy(A);
212:   }
213:   MatDestroy(J);
214:   if (matfdcoloring) {
215:     MatFDColoringDestroy(matfdcoloring);
216:   }
217:   VecDestroy(x);
218:   VecDestroy(r);
219:   SNESDestroy(snes);
220:   DADestroy(user.da);
221:   PetscFinalize();

223:   return(0);
224: }
225: /* ------------------------------------------------------------------- */
228: /* 
229:    FormInitialGuess - Forms initial approximation.

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

235:    Output Parameter:
236:    X - vector
237:  */
238: int FormInitialGuess(AppCtx *user,Vec X)
239: {
240:   int         i,j,Mx,My,ierr,xs,ys,xm,ym;
241:   PetscReal   lambda,temp1,temp,hx,hy;
242:   PetscScalar **x;

245:   DAGetInfo(user->da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
246:                    PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

248:   lambda = user->param;
249:   hx     = 1.0/(PetscReal)(Mx-1);
250:   hy     = 1.0/(PetscReal)(My-1);
251:   temp1  = lambda/(lambda + 1.0);

253:   /*
254:      Get a pointer to vector data.
255:        - For default PETSc vectors, VecGetArray() returns a pointer to
256:          the data array.  Otherwise, the routine is implementation dependent.
257:        - You MUST call VecRestoreArray() when you no longer need access to
258:          the array.
259:   */
260:   DAVecGetArray(user->da,X,&x);

262:   /*
263:      Get local grid boundaries (for 2-dimensional DA):
264:        xs, ys   - starting grid indices (no ghost points)
265:        xm, ym   - widths of local grid (no ghost points)

267:   */
268:   DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);

270:   /*
271:      Compute initial guess over the locally owned part of the grid
272:   */
273:   for (j=ys; j<ys+ym; j++) {
274:     temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
275:     for (i=xs; i<xs+xm; i++) {

277:       if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
278:         /* boundary conditions are all zero Dirichlet */
279:         x[j][i] = 0.0;
280:       } else {
281:         x[j][i] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
282:       }
283:     }
284:   }

286:   /*
287:      Restore vector
288:   */
289:   DAVecRestoreArray(user->da,X,&x);

291:   return(0);
292: }
293: /* ------------------------------------------------------------------- */
296: /* 
297:    FormFunctionLocal - Evaluates nonlinear function, F(x).

299:        Process adiC(36): FormFunctionLocal

301:  */
302: int FormFunctionLocal(DALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
303: {
304:   int         ierr,i,j;
305:   PetscReal   two = 2.0,lambda,hx,hy,hxdhy,hydhx,sc;
306:   PetscScalar u,uxx,uyy;


310:   lambda = user->param;
311:   hx     = 1.0/(PetscReal)(info->mx-1);
312:   hy     = 1.0/(PetscReal)(info->my-1);
313:   sc     = hx*hy*lambda;
314:   hxdhy  = hx/hy;
315:   hydhx  = hy/hx;
316:   /*
317:      Compute function over the locally owned part of the grid
318:   */
319:   for (j=info->ys; j<info->ys+info->ym; j++) {
320:     for (i=info->xs; i<info->xs+info->xm; i++) {
321:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
322:         f[j][i] = x[j][i];
323:       } else {
324:         u       = x[j][i];
325:         uxx     = (two*u - x[j][i-1] - x[j][i+1])*hydhx;
326:         uyy     = (two*u - x[j-1][i] - x[j+1][i])*hxdhy;
327:         f[j][i] = uxx + uyy - sc*PetscExpScalar(u);
328:       }
329:     }
330:   }

332:   PetscLogFlops(11*info->ym*info->xm);
333:   return(0);
334: }

338: /*
339:    FormJacobianLocal - Evaluates Jacobian matrix.


342: */
343: int FormJacobianLocal(DALocalInfo *info,PetscScalar **x,Mat jac,AppCtx *user)
344: {
345:   int          ierr,i,j;
346:   MatStencil   col[5],row;
347:   PetscScalar  lambda,v[5],hx,hy,hxdhy,hydhx,sc;

350:   lambda = user->param;
351:   hx     = 1.0/(PetscReal)(info->mx-1);
352:   hy     = 1.0/(PetscReal)(info->my-1);
353:   sc     = hx*hy*lambda;
354:   hxdhy  = hx/hy;
355:   hydhx  = hy/hx;


358:   /* 
359:      Compute entries for the locally owned part of the Jacobian.
360:       - Currently, all PETSc parallel matrix formats are partitioned by
361:         contiguous chunks of rows across the processors. 
362:       - Each processor needs to insert only elements that it owns
363:         locally (but any non-local elements will be sent to the
364:         appropriate processor during matrix assembly). 
365:       - Here, we set all entries for a particular row at once.
366:       - We can set matrix entries either using either
367:         MatSetValuesLocal() or MatSetValues(), as discussed above.
368:   */
369:   for (j=info->ys; j<info->ys+info->ym; j++) {
370:     for (i=info->xs; i<info->xs+info->xm; i++) {
371:       row.j = j; row.i = i;
372:       /* boundary points */
373:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
374:         v[0] = 1.0;
375:         MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
376:       } else {
377:       /* interior grid points */
378:         v[0] = -hxdhy;                                           col[0].j = j - 1; col[0].i = i;
379:         v[1] = -hydhx;                                           col[1].j = j;     col[1].i = i-1;
380:         v[2] = 2.0*(hydhx + hxdhy) - sc*PetscExpScalar(x[j][i]); col[2].j = row.j; col[2].i = row.i;
381:         v[3] = -hydhx;                                           col[3].j = j;     col[3].i = i+1;
382:         v[4] = -hxdhy;                                           col[4].j = j + 1; col[4].i = i;
383:         MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
384:       }
385:     }
386:   }

388:   /* 
389:      Assemble matrix, using the 2-step process:
390:        MatAssemblyBegin(), MatAssemblyEnd().
391:   */
392:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
393:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
394:   /*
395:      Tell the matrix we will never add a new nonzero location to the
396:      matrix. If we do, it will generate an error.
397:   */
398:   MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR);
399:   return(0);
400: }

402: /*
403:       Variant of FormFunction() that computes the function in Matlab
404: */
405: #if defined(PETSC_HAVE_MATLAB) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE)
406: int FormFunctionMatlab(SNES snes,Vec X,Vec F,void *ptr)
407: {
408:   AppCtx    *user = (AppCtx*)ptr;
409:   int       ierr,Mx,My;
410:   PetscReal lambda,hx,hy;
411:   Vec       localX,localF;
412:   MPI_Comm  comm;

415:   DAGetLocalVector(user->da,&localX);
416:   DAGetLocalVector(user->da,&localF);
417:   PetscObjectSetName((PetscObject)localX,"localX");
418:   PetscObjectSetName((PetscObject)localF,"localF");
419:   DAGetInfo(user->da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
420:                    PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

422:   lambda = user->param;
423:   hx     = 1.0/(PetscReal)(Mx-1);
424:   hy     = 1.0/(PetscReal)(My-1);

426:   PetscObjectGetComm((PetscObject)snes,&comm);
427:   /*
428:      Scatter ghost points to local vector,using the 2-step process
429:         DAGlobalToLocalBegin(),DAGlobalToLocalEnd().
430:      By placing code between these two statements, computations can be
431:      done while messages are in transition.
432:   */
433:   DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
434:   DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);
435:   PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localX);
436:   PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(comm),"localF=ex5m(localX,%18.16e,%18.16e,%18.16e)",hx,hy,lambda);
437:   PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localF);

439:   /*
440:      Insert values into global vector
441:   */
442:   DALocalToGlobal(user->da,localF,INSERT_VALUES,F);
443:   DARestoreLocalVector(user->da,&localX);
444:   DARestoreLocalVector(user->da,&localF);
445:   return(0);
446: }
447: #endif