Actual source code: ex8.c

  1: /* Program usage:  mpiexec -n <procs> ex5 [-help] [all PETSc options] */

  3: static char help[] = "Nonlinear PDE in 2d.\n\
  4: We solve the Bratu equation in a 2D rectangular\n\
  5: domain, using distributed arrays (DMDAs) to partition the parallel grid.\n\n";

  7: /*T
  8:    Concepts: SNES^parallel Bratu example
  9:    Concepts: DMDA^using distributed arrays;
 10:    Processors: n
 11: T*/

 13: /* ------------------------------------------------------------------------

 15:     The Bratu equation is given by the partial differential equation
 16:   
 17:             -alpha*Laplacian u + lambda*e^u = f,  0 < x,y < 1,
 18:   
 19:     with boundary conditions
 20:    
 21:              u = 0  for  x = 0, x = 1, y = 0, y = 1.
 22:   
 23:     A linear finite element approximation is used to discretize the boundary
 24:     value problem on the two triangles which make up each rectangle in the DMDA
 25:     to obtain a nonlinear system of equations.

 27:   ------------------------------------------------------------------------- */

 29: /* 
 30:    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
 31:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
 32:    file automatically includes:
 33:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 34:      petscmat.h - matrices
 35:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 36:      petscviewer.h - viewers               petscpc.h  - preconditioners
 37:      petscksp.h   - linear solvers
 38: */
 39: #include <petscsys.h>
 40: #include <petscbag.h>
 41: #include <petscdmda.h>
 42: #include <petscdmmg.h>
 43: #include <petscsnes.h>

 45: /* 
 46:    User-defined application context - contains data needed by the 
 47:    application-provided call-back routines, FormJacobianLocal() and
 48:    FormFunctionLocal().
 49: */
 50: typedef struct {
 51:    PetscReal alpha;          /* parameter controlling linearity */
 52:    PetscReal lambda;         /* parameter controlling nonlinearity */
 53: } AppCtx;

 55: static PetscScalar Kref[36] = { 0.5,  0.5, -0.5,  0,  0, -0.5,
 56:                                 0.5,  0.5, -0.5,  0,  0, -0.5,
 57:                                -0.5, -0.5,  0.5,  0,  0,  0.5,
 58:                                   0,    0,    0,  0,  0,    0,
 59:                                   0,    0,    0,  0,  0,    0,
 60:                                -0.5, -0.5,  0.5,  0,  0,  0.5};

 62: /* These are */
 63: static PetscScalar quadPoints[8] = {0.17855873, 0.15505103,
 64:                                     0.07503111, 0.64494897,
 65:                                     0.66639025, 0.15505103,
 66:                                     0.28001992, 0.64494897};
 67: static PetscScalar quadWeights[4] = {0.15902069,  0.09097931,  0.15902069,  0.09097931};

 69: /* 
 70:    User-defined routines
 71: */

 80: int main(int argc,char **argv)
 81: {
 82:   DMMG                  *dmmg;                 /* hierarchy manager */
 83:   DM                     da;
 84:   SNES                   snes;                 /* nonlinear solver */
 85:   AppCtx                *user;                 /* user-defined work context */
 86:   PetscBag               bag;
 87:   PetscInt               its;                  /* iterations for convergence */
 88:   SNESConvergedReason    reason;
 89:   PetscBool              drawContours;         /* flag for drawing contours */
 90:   PetscErrorCode         ierr;
 91:   PetscReal              lambda_max = 6.81, lambda_min = 0.0, error;

 93:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 94:      Initialize program
 95:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 96:   PetscInitialize(&argc,&argv,(char *)0,help);

 98:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 99:      Initialize problem parameters
100:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
101:   PetscBagCreate(PETSC_COMM_WORLD, sizeof(AppCtx), &bag);
102:   PetscBagGetData(bag, (void **) &user);
103:   PetscBagSetName(bag, "params", "Parameters for SNES example 4");
104:   PetscBagRegisterReal(bag, &user->alpha, 1.0, "alpha", "Linear coefficient");
105:   PetscBagRegisterReal(bag, &user->lambda, 6.0, "lambda", "Nonlinear coefficient");
106:   PetscBagSetFromOptions(bag);
107:   PetscOptionsGetReal(PETSC_NULL,"-alpha",&user->alpha,PETSC_NULL);
108:   PetscOptionsGetReal(PETSC_NULL,"-lambda",&user->lambda,PETSC_NULL);
109:   if (user->lambda > lambda_max || user->lambda < lambda_min) {
110:     SETERRQ3(PETSC_COMM_SELF,1,"Lambda %G is out of range [%G, %G]", user->lambda, lambda_min, lambda_max);
111:   }

113:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114:      Create multilevel DM data structure (DMMG) to manage hierarchical solvers
115:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
116:   DMMGCreate(PETSC_COMM_WORLD,1,user,&dmmg);

118:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119:      Create distributed array (DA) to manage parallel grid and vectors
120:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
121:   DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,-3,-3,PETSC_DECIDE,PETSC_DECIDE,
122:                     1,1,PETSC_NULL,PETSC_NULL,&da);
123:   DMDASetFieldName(da, 0, "ooblek");
124:   DMMGSetDM(dmmg, (DM) da);
125:   DMDestroy(&da);

127:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128:      Set the discretization functions
129:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130:   DMMGSetSNESLocal(dmmg, FormFunctionLocal, FormJacobianLocal, 0, 0);
131:   DMMGSetFromOptions(dmmg);

133:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134:      Evaluate initial guess
135:      Note: The user should initialize the vector, x, with the initial guess
136:      for the nonlinear solver prior to calling SNESSolve().  In particular,
137:      to employ an initial guess of zero, the user should explicitly set
138:      this vector to zero by calling VecSet().
139:   */
140:   DMMGSetInitialGuess(dmmg, FormInitialGuess);

142:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143:      Solve nonlinear system
144:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
145:   DMMGSolve(dmmg);
146:   snes = DMMGGetSNES(dmmg);
147:   SNESGetIterationNumber(snes,&its);
148:   SNESGetConvergedReason(snes, &reason);

150:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
151:   PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %D, %s\n",its,SNESConvergedReasons[reason]);
152:   L_2Error(DMMGGetDM(dmmg), DMMGGetx(dmmg), &error, user);
153:   PetscPrintf(PETSC_COMM_WORLD,"L_2 error in the solution: %G\n", error);
154:   PrintVector(dmmg[0], DMMGGetx(dmmg));

156:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157:      Free work space.  All PETSc objects should be destroyed when they
158:      are no longer needed.
159:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
160:   DMMGDestroy(dmmg);
161:   PetscBagDestroy(&bag);
162:   PetscFinalize();
163:   return(0);
164: }

168: PetscErrorCode PrintVector(DMMG dmmg, Vec U)
169: {
170:   DM             da =  dmmg->dm;
171:   PetscScalar  **u;
172:   PetscInt       i,j,xs,ys,xm,ym;

176:   DMDAVecGetArray(da,U,&u);
177:   DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
178:   for(j = ys+ym-1; j >= ys; j--) {
179:     for(i = xs; i < xs+xm; i++) {
180:       printf("u[%d][%d] = %G ", j, i, u[j][i]);
181:     }
182:     printf("\n");
183:   }
184:   DMDAVecRestoreArray(da,U,&u);
185:   return(0);
186: }

190: PetscErrorCode ExactSolution(PetscReal x, PetscReal y, PetscScalar *u)
191: {
193:   *u = x*x;
194:   return(0);
195: }

199: /* 
200:    FormInitialGuess - Forms initial approximation.

202:    Input Parameters:
203:    dmmg - The DMMG context
204:    X - vector

206:    Output Parameter:
207:    X - vector
208: */
209: PetscErrorCode FormInitialGuess(DMMG dmmg,Vec X)
210: {
211:   AppCtx        *user = (AppCtx *) dmmg->user;
212:   DM             da =  dmmg->dm;
213:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
215:   PetscReal      lambda,temp1,temp,hx,hy;
216:   PetscScalar    **x;

219:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
220:                    PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

222:   lambda = user->lambda;
223:   hx     = 1.0/(PetscReal)(Mx-1);
224:   hy     = 1.0/(PetscReal)(My-1);
225:   if (lambda == 0.0) {
226:     temp1  = 0.0;
227:   } else {
228:     temp1  = lambda/(lambda + 1.0);
229:   }

231:   /*
232:      Get a pointer to vector data.
233:        - For default PETSc vectors, VecGetArray() returns a pointer to
234:          the data array.  Otherwise, the routine is implementation dependent.
235:        - You MUST call VecRestoreArray() when you no longer need access to
236:          the array.
237:   */
238:   DMDAVecGetArray(da,X,&x);

240:   /*
241:      Get local grid boundaries (for 2-dimensional DMDA):
242:        xs, ys   - starting grid indices (no ghost points)
243:        xm, ym   - widths of local grid (no ghost points)

245:   */
246:   DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);

248:   /*
249:      Compute initial guess over the locally owned part of the grid
250:   */
251:   for (j=ys; j<ys+ym; j++) {
252:     temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
253:     for (i=xs; i<xs+xm; i++) {

255:       if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
256:         /* boundary conditions are all zero Dirichlet */
257:         x[j][i] = 0.0;
258:       } else {
259:         x[j][i] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
260:       }
261:     }
262:   }

264:   DMDAVecRestoreArray(da,X,&x);
265:   PrintVector(dmmg, X);
266:   return(0);
267: }

271: PetscErrorCode constantResidual(PetscReal lambda, PetscBool  isLower, int i, int j, PetscReal hx, PetscReal hy, PetscScalar r[])
272: {
273:   PetscScalar rLocal[3] = {0.0, 0.0, 0.0};
274:   PetscScalar phi[3] = {0.0, 0.0, 0.0};
275:   PetscReal   xI = i*hx, yI = j*hy, hxhy = hx*hy, x, y;
276:   PetscScalar res;
277:   PetscInt    q, k;

280:   for(q = 0; q < 4; q++) {
281:     phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
282:     phi[1] = quadPoints[q*2];
283:     phi[2] = quadPoints[q*2+1];
284:     /* These are currently wrong */
285:     x      = xI + quadPoints[q*2]*hx;
286:     y      = yI + quadPoints[q*2+1]*hy;
287:     res    = quadWeights[q]*(2.0);
288:     for(k = 0; k < 3; k++) {
289:       rLocal[k] += phi[k]*res;
290:     }
291:   }
292:   for(k = 0; k < 3; k++) {
293:     printf("  constLocal[%d] = %g\n", k, lambda*hxhy*rLocal[k]);
294:     r[k] += lambda*hxhy*rLocal[k];
295:   }
296:   return(0);
297: }

299: PetscErrorCode nonlinearResidual(PetscReal lambda, PetscScalar u[], PetscScalar r[]) {
300:   PetscScalar rLocal[3] = {0.0, 0.0, 0.0};
301:   PetscScalar phi[3] = {0.0, 0.0, 0.0};
302:   PetscScalar res;
303:   PetscInt q;

306:   for(q = 0; q < 4; q++) {
307:     phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
308:     phi[1] = quadPoints[q*2];
309:     phi[2] = quadPoints[q*2+1];
310:     res    = quadWeights[q]*PetscExpScalar(u[0]*phi[0]+ u[1]*phi[1] + u[2]*phi[2]+ u[3]*phi[3]);
311:     rLocal[0] += phi[0]*res;
312:     rLocal[1] += phi[1]*res;
313:     rLocal[2] += phi[2]*res;
314:   }
315:   r[0] += lambda*rLocal[0];
316:   r[1] += lambda*rLocal[1];
317:   r[2] += lambda*rLocal[2];
318:   return(0);
319: }

323: /* 
324:    FormFunctionLocal - Evaluates nonlinear function, F(x).

326:        Process adiC(36): FormFunctionLocal

328:  */
329: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
330: {
331:   PetscScalar    uLocal[3];
332:   PetscScalar    rLocal[3];
333:   PetscScalar    G[4];
334:   PetscScalar    uExact;
335:   PetscReal      alpha,lambda,hx,hy,hxhy,sc,detJInv;
336:   PetscInt       i,j,k,l;


341:   /* Naive Jacobian calculation:

343:      J = / 1/hx  0   \ J^{-1} = / hx   0 \  1/|J| = hx*hy = |J^{-1}|
344:          \  0   1/hy /          \  0  hy /
345:    */
346:   alpha   = user->alpha;
347:   lambda  = user->lambda;
348:   hx      = 1.0/(PetscReal)(info->mx-1);
349:   hy      = 1.0/(PetscReal)(info->my-1);
350:   sc      = hx*hy*lambda;
351:   hxhy    = hx*hy;
352:   detJInv = hxhy;
353:   G[0] = (1.0/(hx*hx)) * detJInv;
354:   G[1] = 0.0;
355:   G[2] = G[1];
356:   G[3] = (1.0/(hy*hy)) * detJInv;
357:   for(k = 0; k < 4; k++) {
358:     printf("G[%d] = %g\n", k, G[k]);
359:   }

361:   /* Zero the vector */
362:   PetscMemzero((void *) &(f[info->xs][info->ys]), info->xm*info->ym*sizeof(PetscScalar));
363:   /* Compute function over the locally owned part of the grid. For each
364:      vertex (i,j), we consider the element below:

366:        2 (1)    (0)
367:      i,j+1 --- i+1,j+1
368:        |  \      |
369:        |   \     |
370:        |    \    |
371:        |     \   |
372:        |      \  |
373:       i,j  --- i+1,j
374:        0         1 (2)

376:      and therefore we do not loop over the last vertex in each dimension.
377:   */
378:   for(j = info->ys; j < info->ys+info->ym-1; j++) {
379:     for(i = info->xs; i < info->xs+info->xm-1; i++) {
380:       /* Lower element */
381:       uLocal[0] = x[j][i];
382:       uLocal[1] = x[j][i+1];
383:       uLocal[2] = x[j+1][i];
384:       printf("Solution ElementVector for (%d, %d)\n", i, j);
385:       for(k = 0; k < 3; k++) {
386:         printf("  uLocal[%d] = %g\n", k, uLocal[k]);
387:       }
388:       for(k = 0; k < 3; k++) {
389:         rLocal[k] = 0.0;
390:         for(l = 0; l < 3; l++) {
391:           rLocal[k] += (G[0]*Kref[(k*2*3 + l)*2]+G[1]*Kref[(k*2*3 + l)*2+1]+G[2]*Kref[((k*2+1)*3 + l)*2]+G[3]*Kref[((k*2+1)*3 + l)*2+1])*uLocal[l];
392:         }
393:         rLocal[k] *= alpha;
394:       }
395:       printf("Laplacian ElementVector for (%d, %d)\n", i, j);
396:       for(k = 0; k < 3; k++) {
397:         printf("  rLocal[%d] = %g\n", k, rLocal[k]);
398:       }
399:       constantResidual(1.0, PETSC_TRUE, i, j, hx, hy, rLocal);
400:       printf("Laplacian+Constant ElementVector for (%d, %d)\n", i, j);
401:       for(k = 0; k < 3; k++) {
402:         printf("  rLocal[%d] = %g\n", k, rLocal[k]);
403:       }
404:       nonlinearResidual(0.0*sc, uLocal, rLocal);
405:       printf("Full nonlinear ElementVector for (%d, %d)\n", i, j);
406:       for(k = 0; k < 3; k++) {
407:         printf("  rLocal[%d] = %g\n", k, rLocal[k]);
408:       }
409:       f[j][i]   += rLocal[0];
410:       f[j][i+1] += rLocal[1];
411:       f[j+1][i] += rLocal[2];
412:       /* Upper element */
413:       uLocal[0] = x[j+1][i+1];
414:       uLocal[1] = x[j+1][i];
415:       uLocal[2] = x[j][i+1];
416:       printf("Solution ElementVector for (%d, %d)\n", i, j);
417:       for(k = 0; k < 3; k++) {
418:         printf("  uLocal[%d] = %g\n", k, uLocal[k]);
419:       }
420:       for(k = 0; k < 3; k++) {
421:         rLocal[k] = 0.0;
422:         for(l = 0; l < 3; l++) {
423:           rLocal[k] += (G[0]*Kref[(k*2*3 + l)*2]+G[1]*Kref[(k*2*3 + l)*2+1]+G[2]*Kref[((k*2+1)*3 + l)*2]+G[3]*Kref[((k*2+1)*3 + l)*2+1])*uLocal[l];
424:         }
425:         rLocal[k] *= alpha;
426:       }
427:       printf("Laplacian ElementVector for (%d, %d)\n", i, j);
428:       for(k = 0; k < 3; k++) {
429:         printf("  rLocal[%d] = %g\n", k, rLocal[k]);
430:       }
431:       constantResidual(1.0, PETSC_BOOL, i, j, hx, hy, rLocal);
432:       printf("Laplacian+Constant ElementVector for (%d, %d)\n", i, j);
433:       for(k = 0; k < 3; k++) {
434:         printf("  rLocal[%d] = %g\n", k, rLocal[k]);
435:       }
436:       nonlinearResidual(0.0*sc, uLocal, rLocal);
437:       printf("Full nonlinear ElementVector for (%d, %d)\n", i, j);
438:       for(k = 0; k < 3; k++) {
439:         printf("  rLocal[%d] = %g\n", k, rLocal[k]);
440:       }
441:       f[j+1][i+1] += rLocal[0];
442:       f[j+1][i]   += rLocal[1];
443:       f[j][i+1]   += rLocal[2];
444:       /* Boundary conditions */
445:       if (i == 0 || j == 0) {
446:         ExactSolution(i*hx, j*hy, &uExact);
447:         f[j][i] = x[j][i] - uExact;
448:       }
449:       if ((i == info->mx-2) || (j == 0)) {
450:         ExactSolution((i+1)*hx, j*hy, &uExact);
451:         f[j][i+1] = x[j][i+1] - uExact;
452:       }
453:       if ((i == info->mx-2) || (j == info->my-2)) {
454:         ExactSolution((i+1)*hx, (j+1)*hy, &uExact);
455:         f[j+1][i+1] = x[j+1][i+1] - uExact;
456:       }
457:       if ((i == 0) || (j == info->my-2)) {
458:         ExactSolution(i*hx, (j+1)*hy, &uExact);
459:         f[j+1][i] = x[j+1][i] - uExact;
460:       }
461:     }
462:   }

464:   for(j = info->ys+info->ym-1; j >= info->ys; j--) {
465:     for(i = info->xs; i < info->xs+info->xm; i++) {
466:       printf("f[%d][%d] = %g ", j, i, f[j][i]);
467:     }
468:     printf("\n");
469:   }
470:   PetscLogFlops(68.0*(info->ym-1)*(info->xm-1));
471:   return(0);
472: }

474: PetscErrorCode nonlinearJacobian(PetscReal lambda, PetscScalar u[], PetscScalar J[]) {
476:   return(0);
477: }

481: /*
482:    FormJacobianLocal - Evaluates Jacobian matrix.
483: */
484: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info,PetscScalar **x,Mat jac,AppCtx *user)
485: {
486:   PetscScalar    JLocal[16], uLocal[4];
487:   MatStencil     rows[4], cols[4], ident;
488:   PetscInt       lowerRow[3] = {0, 1, 3};
489:   PetscInt       upperRow[3] = {2, 3, 1};
490:   PetscInt       hasLower[3], hasUpper[3], localRows[4];
491:   PetscScalar    alpha,lambda,hx,hy,hxhy,detJInv,G[4],sc,one = 1.0;
492:   PetscInt       i,j,k,l,numRows;

496:   alpha  = user->alpha;
497:   lambda = user->lambda;
498:   hx     = 1.0/(PetscReal)(info->mx-1);
499:   hy     = 1.0/(PetscReal)(info->my-1);
500:   sc     = hx*hy*lambda;
501:   hxhy   = hx*hy;
502:   detJInv = hxhy;
503:   G[0] = (1.0/(hx*hx)) * detJInv;
504:   G[1] = 0.0;
505:   G[2] = G[1];
506:   G[3] = (1.0/(hy*hy)) * detJInv;
507:   for(k = 0; k < 4; k++) {
508:     printf("G[%d] = %g\n", k, G[k]);
509:   }

511:   MatZeroEntries(jac);
512:   /* 
513:      Compute entries for the locally owned part of the Jacobian.
514:       - Currently, all PETSc parallel matrix formats are partitioned by
515:         contiguous chunks of rows across the processors. 
516:       - Each processor needs to insert only elements that it owns
517:         locally (but any non-local elements will be sent to the
518:         appropriate processor during matrix assembly). 
519:       - Here, we set all entries for a particular row at once.
520:       - We can set matrix entries either using either
521:         MatSetValuesLocal() or MatSetValues(), as discussed above.
522:   */
523:   for (j=info->ys; j<info->ys+info->ym-1; j++) {
524:     for (i=info->xs; i<info->xs+info->xm-1; i++) {
525:       PetscMemzero(JLocal, 16 * sizeof(PetscScalar));
526:       numRows = 0;
527:       /* Lower element */
528:       uLocal[0] = x[j][i];
529:       uLocal[1] = x[j][i+1];
530:       uLocal[2] = x[j+1][i+1];
531:       uLocal[3] = x[j+1][i];
532:       /* i,j */
533:       if (i == 0 || j == 0) {
534:         hasLower[0] = 0;
535:         ident.i = i; ident.j = j;
536:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
537:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
538:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
539:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
540:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
541:       } else {
542:         hasLower[0] = 1;
543:         localRows[0] = numRows;
544:         rows[numRows].i = i; rows[numRows].j = j;
545:         numRows++;
546:       }
547:       cols[0].i = i; cols[0].j = j;
548:       /* i+1,j */
549:       if ((i == info->mx-2) || (j == 0)) {
550:         hasLower[1] = 0;
551:         hasUpper[2] = 0;
552:         ident.i = i+1; ident.j = j;
553:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
554:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
555:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
556:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
557:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
558:       } else {
559:         localRows[1] = numRows;
560:         hasLower[1] = 1;
561:         hasUpper[2] = 1;
562:         localRows[1] = numRows;
563:         rows[numRows].i = i+1; rows[numRows].j = j;
564:         numRows++;
565:       }
566:       cols[1].i = i+1; cols[1].j = j;
567:       /* i+1,j+1 */
568:       if ((i == info->mx-2) || (j == info->my-2)) {
569:         hasUpper[0] = 0;
570:         ident.i = i+1; ident.j = j+1;
571:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
572:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
573:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
574:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
575:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
576:       } else {
577:         hasUpper[0] = 1;
578:         localRows[2] = numRows;
579:         rows[numRows].i = i+1; rows[numRows].j = j+1;
580:         numRows++;
581:       }
582:       cols[2].i = i+1; cols[2].j = j+1;
583:       /* i,j+1 */
584:       if ((i == 0) || (j == info->my-2)) {
585:         hasLower[2] = 0;
586:         hasUpper[1] = 0;
587:         ident.i = i; ident.j = j+1;
588:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
589:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
590:         MatSetValuesStencil(jac,1,&ident,1,&ident,&one,INSERT_VALUES);
591:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
592:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
593:       } else {
594:         hasLower[2] = 1;
595:         hasUpper[1] = 1;
596:         localRows[3] = numRows;
597:         rows[numRows].i = i; rows[numRows].j = j+1;
598:         numRows++;
599:       }
600:       cols[3].i = i; cols[3].j = j+1;

602:       /* Lower Element */
603:       for(k = 0; k < 3; k++) {
604:         if (!hasLower[k]) continue;
605:         for(l = 0; l < 3; l++) {
606:           JLocal[localRows[lowerRow[k]]*4 + lowerRow[l]] += alpha*(G[0]*Kref[(k*2*3 + l)*2]+G[1]*Kref[(k*2*3 + l)*2+1]+G[2]*Kref[((k*2+1)*3 + l)*2]+G[3]*Kref[((k*2+1)*3 + l)*2+1]);
607:         }
608:       }
609:       /* Upper Element */
610:       for(k = 0; k < 3; k++) {
611:         if (!hasUpper[k]) continue;
612:         for(l = 0; l < 3; l++) {
613:           JLocal[localRows[upperRow[k]]*4 + upperRow[l]] += alpha*(G[0]*Kref[(k*2*3 + l)*2]+G[1]*Kref[(k*2*3 + l)*2+1]+G[2]*Kref[((k*2+1)*3 + l)*2]+G[3]*Kref[((k*2+1)*3 + l)*2+1]);
614:         }
615:       }

617:       nonlinearJacobian(-1.0*sc, uLocal, JLocal);
618:       printf("Element matrix for (%d, %d)\n", i, j);
619:       printf("   col  ");
620:       for(l = 0; l < 4; l++) {
621:         printf("(%d, %d) ", cols[l].i, cols[l].j);
622:       }
623:       printf("\n");
624:       for(k = 0; k < numRows; k++) {
625:         printf("row (%d, %d): ", rows[k].i, rows[k].j);
626:         for(l = 0; l < 4; l++) {
627:           printf("%8.6g ", JLocal[k*4 + l]);
628:         }
629:         printf("\n");
630:       }
631:       MatSetValuesStencil(jac,numRows,rows,4,cols,JLocal,ADD_VALUES);
632:     }
633:   }

635:   /* 
636:      Assemble matrix, using the 2-step process:
637:        MatAssemblyBegin(), MatAssemblyEnd().
638:   */
639:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
640:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
641:   /*
642:      Tell the matrix we will never add a new nonzero location to the
643:      matrix. If we do, it will generate an error.
644:   */
645:   MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
646:   return(0);
647: }

651: /*
652:   L_2Error - Integrate the L_2 error of our solution over each face
653: */
654: PetscErrorCode L_2Error(DM da, Vec fVec, double *error, AppCtx *user)
655: {
656:   DMDALocalInfo info;
657:   Vec fLocalVec;
658:   PetscScalar **f;
659:   PetscScalar u, uExact, uLocal[4];
660:   PetscScalar hx, hy, hxhy, x, y, phi[3];
661:   PetscInt i, j, q;

665:   DMDAGetLocalInfo(da, &info);
666:   DMGetLocalVector(da, &fLocalVec);
667:   DMGlobalToLocalBegin(da,fVec, INSERT_VALUES, fLocalVec);
668:   DMGlobalToLocalEnd(da,fVec, INSERT_VALUES, fLocalVec);
669:   DMDAVecGetArray(da, fLocalVec, &f);

671:   *error = 0.0;
672:   hx     = 1.0/(PetscReal)(info.mx-1);
673:   hy     = 1.0/(PetscReal)(info.my-1);
674:   hxhy   = hx*hy;
675:   for(j = info.ys; j < info.ys+info.ym-1; j++) {
676:     for(i = info.xs; i < info.xs+info.xm-1; i++) {
677:       uLocal[0] = f[j][i];
678:       uLocal[1] = f[j][i+1];
679:       uLocal[2] = f[j+1][i+1];
680:       uLocal[3] = f[j+1][i];
681:       /* Lower element */
682:       for(q = 0; q < 4; q++) {
683:         phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
684:         phi[1] = quadPoints[q*2];
685:         phi[2] = quadPoints[q*2+1];
686:         u = uLocal[0]*phi[0] + uLocal[1]*phi[1] + uLocal[3]*phi[2];
687:         x = (quadPoints[q*2] + i)*hx;
688:         y = (quadPoints[q*2+1] + j)*hy;
689:         ExactSolution(x, y, &uExact);
690:         *error += hxhy*quadWeights[q]*((u - uExact)*(u - uExact));
691:       }
692:       /* Upper element */
693:       /*
694:         The affine map from the lower to the upper is

696:         / x_U \ = / -1  0 \ / x_L \ + / hx \
697:         \ y_U /   \  0 -1 / \ y_L /   \ hy /
698:        */
699:       for(q = 0; q < 4; q++) {
700:         phi[0] = 1.0 - quadPoints[q*2] - quadPoints[q*2+1];
701:         phi[1] = quadPoints[q*2];
702:         phi[2] = quadPoints[q*2+1];
703:         u = uLocal[2]*phi[0] + uLocal[3]*phi[1] + uLocal[1]*phi[2];
704:         x = (1.0 - quadPoints[q*2] + i)*hx;
705:         y = (1.0 - quadPoints[q*2+1] + j)*hy;
706:         ExactSolution(x, y, &uExact);
707:         *error += hxhy*quadWeights[q]*((u - uExact)*(u - uExact));
708:       }
709:     }
710:   }

712:   DMDAVecRestoreArray(da, fLocalVec, &f);
713:   /* DMLocalToGlobalBegin(da,xLocalVec,ADD_VALUES,xVec); */
714:   /* DMLocalToGlobalEnd(da,xLocalVec,ADD_VALUES,xVec); */
715:   DMRestoreLocalVector(da, &fLocalVec);
716:   return(0);
717: }