Actual source code: poisson2.c

  1: #ifdef PETSC_RCS_HEADER
  2: static char vcid[] = "$Id: poisson2.c,v 1.5 2000/02/01 17:02:51 knepley Exp $";
  3: #endif

  5: static char help[] = "This is the Poisson problem with Dirchlet boundary conditions.\n\n";

  7: int POISSON_COOKIE;
  8: int POISSON_ComputeField;

 10: #include "poisson2.h"

 12: /*
 13:   Here we are solving the Poisson equation,

 15:     -\Delta u = f

 17:   but we introduce an auxiliary variable:

 19:     v = -\nabla u

 21:   so that the Poisson equation becomes

 23:          -v        -   \nabla u = 0           /    -I         -\nabla \   / 0 \
 24:                                        ===>   |                       | = |   |
 25:     \nabla \cdot v              = f           \ \nabla\cdot      0    /   \ f /

 27:   Thus we now let

 29:     SolutionFunction()         --> u_S (which we prescribe)
 30:     GradientSolutionFunction() --> v_S (which we prescribe)
 31:     RhsFunction()              --> f
 32: */
 33: #undef  __FUNCT__
 35: int main(int argc, char **args) {
 36:   PoissonContext ctx; /* Holds problem specific information */
 37:   int           ierr;

 40:   PetscInitialize(&argc, &args, 0, help);                                    CHKERRABORT(PETSC_COMM_WORLD, ierr);

 42:   PoissonContextCreate(PETSC_COMM_WORLD, &ctx);                              CHKERRABORT(PETSC_COMM_WORLD, ierr);
 43:   PoissonComputeField(ctx);                                                  CHKERRABORT(PETSC_COMM_WORLD, ierr);
 44:   PoissonContextDestroy(ctx);                                                CHKERRABORT(PETSC_COMM_WORLD, ierr);

 46:   CHKMEMQ;
 47:   PetscFinalize();
 48:   return(0);
 49: }

 53: /*@C
 54:   PoissonInitializePackage - This function initializes everything in the Poisson package.

 56:   Input Parameter:
 57:   path - The dynamic library path, or PETSC_NULL

 59:   Level: developer

 61: .keywords: Poisson, initialize, package
 62: .seealso: PetscInitialize()
 63: @*/
 64: int PoissonInitializePackage(char *path) {
 65:   static PetscTruth initialized = PETSC_FALSE;
 66:   char              logList[256];
 67:   char             *className;
 68:   PetscTruth        opt;
 69:   int               ierr;

 72:   if (initialized == PETSC_TRUE) return(0);
 73:   initialized = PETSC_TRUE;
 74:   /* Register Classes */
 75:   PetscLogClassRegister(&POISSON_COOKIE, "Poisson");
 76:   /* Register Constructors and Serializers */
 77:   /* Register Events */
 78:   PetscLogEventRegister(&POISSON_ComputeField, "PoissonComputeField", POISSON_COOKIE);
 79:   /* Process info exclusions */
 80:   PetscOptionsGetString(PETSC_NULL, "-log_info_exclude", logList, 256, &opt);
 81:   if (opt == PETSC_TRUE) {
 82:     PetscStrstr(logList, "PoissonContext", &className);
 83:     if (className) {
 84:       PetscLogInfoDeactivateClass(POISSON_COOKIE);
 85:     }
 86:   }
 87:   /* Process summary exclusions */
 88:   PetscOptionsGetString(PETSC_NULL, "-log_summary_exclude", logList, 256, &opt);
 89:   if (opt == PETSC_TRUE) {
 90:     PetscStrstr(logList, "PoissonContext", &className);
 91:     if (className) {
 92:       PetscLogEventDeactivateClass(POISSON_COOKIE);
 93:     }
 94:   }
 95:   return(0);
 96: }

 98: /*-------------------------------------------- PoissonContext Creation ------------------------------------------------*/
 99: #undef  __FUNCT__
101: /*@
102:   PoissonContextCreate - This function initializes the Poisson context.

104:   Collective on MPI_Comm

106:   Input Parameter:
107: . comm - The communicator

109:   Output Parameter:
110: . sCtx - The PoissonContext

112:   Level: beginner

114: .keywords: Poisson, context, create
115: .seealso: PoissonContextDestroy(), PoissonContextPrint(), PoissonContextSetup()
116: @*/
117: int PoissonContextCreate(MPI_Comm comm, PoissonContext *sCtx) {
118:   PoissonContext ctx;

121:   /* Setup context */
122:   PetscHeaderCreate(ctx, _PoissonContext, int, PETSC_VIEWER_COOKIE, -1, "PoissonContext", comm, 0, 0);
123:   PetscLogObjectCreate(ctx);
124:   PetscLogObjectMemory(ctx, sizeof(struct _PoissonContext));

126:   /* Initialize subobjects */
127:   ctx->grid       = PETSC_NULL;
128:   ctx->sles       = PETSC_NULL;
129:   ctx->A          = PETSC_NULL;
130:   ctx->u          = PETSC_NULL;
131:   ctx->f          = PETSC_NULL;
132:   ctx->uExact     = PETSC_NULL;
133:   /* Setup domain */
134:   ctx->geometryCtx.size[0]  = 2.0;
135:   ctx->geometryCtx.size[1]  = 2.0;
136:   ctx->geometryCtx.start[0] = 0.0;
137:   ctx->geometryCtx.start[1] = 0.0;
138:   /* Setup refinement */
139:   ctx->geometryCtx.maxArea  = 0.5;
140:   ctx->geometryCtx.areaFunc = PointFunctionConstant;
141:   ctx->geometryCtx.areaCtx  = PETSC_NULL;
142:   /* Initialize problem loop */
143:   ctx->dim        = 2;
144:   ctx->linear     = PETSC_FALSE;
145:   ctx->flipDisc   = PETSC_FALSE;
146:   ctx->refineStep = 0;
147:   ctx->numLoops   = 0;

149:   *sCtx = ctx;
150:   return(0);
151: }

153: #undef  __FUNCT__
155: /*@
156:   PoissonContextDestroy - This function destroys the Poisson context.

158:   Collective on PoissonContext

160:   Input Parameter:
161: . ctx - The PoissonContext

163:   Level: beginner

165: .keywords: Poisson, context, destroy
166: .seealso: PoissonContextCreate(), PoissonContextSetup()
167: @*/
168: int PoissonContextDestroy(PoissonContext ctx) {
169:   Grid grid = ctx->grid;
170:   int  ierr;

173:   if (--ctx->refct > 0) SETERRQ(PETSC_ERR_PLIB, "Poisson context should not be referenced more than once");
174:   GridFinalizeBC(grid);
175:   GridDestroy(grid);
176:   PetscLogObjectDestroy(ctx);
177:   PetscHeaderDestroy(ctx);
178:   return(0);
179: }

181: /*--------------------------------------------- PoissonContext Setup -------------------------------------------------*/
182: #undef  __FUNCT__
184: int PoissonContextSetup(PoissonContext ctx) {
185:   PetscTruth opt;
186:   int        ierr;

189:   /* Determine the problem dimension */
190:   PetscOptionsGetInt(PETSC_NULL, "-dim", &ctx->dim, &opt);
191:   /* Determine the element type */
192:   PetscOptionsHasName(PETSC_NULL, "-linear", &ctx->linear);
193:   /* Determine the discretization smoothness */
194:   PetscOptionsHasName(PETSC_NULL, "-flip_disc", &ctx->flipDisc);
195:   /* Determine how many systems to solve */
196:   PetscOptionsGetInt(PETSC_NULL, "-num_systems", &ctx->numLoops, &opt);
197:   /* The first iteration at which to refine the mesh */
198:   PetscOptionsGetInt(PETSC_NULL, "-refine_step", &ctx->refineStep, &opt);
199:   /* The maximum area of any triangle in the refined mesh */
200:   PetscOptionsGetReal("mesh", "-max_area", &ctx->geometryCtx.maxArea, &opt);

202:   /* Create main problem */
203:   PoissonContextCreateGrid(ctx);
204:   return(0);
205: }

207: /*------------------------------------------------ Grid Creation -----------------------------------------------------*/
208: #undef  __FUNCT__
210: /*@
211:   PoissonContextCreateMeshBoundary - This function creates a mesh boundary for the main problem.

213:   Collective on PoissonContext

215:   Input Parameter:
216: . ctx - A PoissonContext with problem specific information

218:   Level: beginner

220: .seealso PoissonContextDestroyMeshBoundary(), PoissonContextCreateMesh()
221: @*/
222: int PoissonContextCreateMeshBoundary(PoissonContext ctx) {
223:   MPI_Comm             comm;
224:   MeshGeometryContext *geomCtx = &ctx->geometryCtx;
225:   int                  ierr;

228:   PetscObjectGetComm((PetscObject) ctx, &comm);
229:   switch(ctx->dim) {
230:   case 1:
231:     MeshBoundary1DCreateSimple(comm, geomCtx, &ctx->boundaryCtx);
232:     break;
233:   case 2:
234:     MeshBoundary2DCreateSimple(comm, geomCtx, &ctx->boundaryCtx);
235:     break;
236:   default:
237:     SETERRQ1(PETSC_ERR_SUP, "Do not support meshes of dimension %d", ctx->dim);
238:   }
239:   return(0);
240: }

242: #undef  __FUNCT__
244: /*@
245:   PoissonContextDestroyMeshBoundary - This function destroys the mesh boundary for the main problem.

247:   Collective on PoissonContext

249:   Input Parameter:
250: . ctx - A PoissonContext with problem specific information

252:   Level: beginner

254: .seealso PoissonContextCreateMeshBoundary(), PoissonContextCreateMesh()
255: @*/
256: int PoissonContextDestroyMeshBoundary(PoissonContext ctx) {
257:   MPI_Comm comm;
258:   int      ierr;

261:   PetscObjectGetComm((PetscObject) ctx, &comm);
262:   switch(ctx->dim) {
263:   case 1:
264:     MeshBoundary1DDestroy(comm, &ctx->boundaryCtx);
265:     break;
266:   case 2:
267:     MeshBoundary2DDestroy(comm, &ctx->boundaryCtx);
268:     break;
269:   default:
270:     SETERRQ1(PETSC_ERR_SUP, "Do not support meshes of dimension %d", ctx->dim);
271:   }
272:   return(0);
273: }

275: #undef  __FUNCT__
277: /*@
278:   PoissonContextCreateMesh - This function creates a mesh for the main problem.

280:   Collective on PoissonContext

282:   Input Parameter:
283: . ctx - A PoissonContext with problem specific information

285:   Output Parameter:
286: . m   - The Mesh

288:   Options Database Keys:
289: . mesh_refine   - Refines the mesh based on area criteria
290: . mesh_max_area - The maximum area of an element

292:   Level: beginner

294: .seealso PoissonRefineGrid(), PoissonDestroyGrid()
295: @*/
296: int PoissonContextCreateMesh(PoissonContext ctx, Mesh *m) {
297:   MeshGeometryContext *geomCtx = &ctx->geometryCtx;
298:   MPI_Comm             comm;
299:   Mesh                 mesh;
300:   Partition            part;
301:   int                  totalElements, totalNodes, totalEdges;
302:   char                 name[1024];
303:   int                  d;
304:   int                  ierr;

307:   PetscObjectGetComm((PetscObject) ctx, &comm);
308:   PoissonContextCreateMeshBoundary(ctx);
309:   MeshCreate(comm, &mesh);
310:   MeshSetDimension(mesh, ctx->dim);
311:   for(d = 0; d < ctx->dim; d++) {
312:     MeshSetPeriodicDimension(mesh, d, geomCtx->isPeriodic[d]);
313:   }
314:   switch(ctx->dim) {
315:   case 1:
316:     MeshSetNumCorners(mesh, 3);
317:     break;
318:   case 2:
319:     if (ctx->linear == PETSC_TRUE) {
320:       MeshSetNumCorners(mesh, 4);
321:     } else {
322:       MeshSetNumCorners(mesh, 6);
323:     }
324:     break;
325:   default:
326:     SETERRQ1(PETSC_ERR_SUP, "Do not support meshes of dimension %d", ctx->dim);
327:   }
328:   MeshSetBoundary(mesh, &ctx->boundaryCtx);
329:   sprintf(name, "mesh.r%.4g", geomCtx->maxArea);
330:   PetscObjectSetName((PetscObject) mesh, name);
331:   MeshSetFromOptions(mesh);
332:   PoissonContextDestroyMeshBoundary(ctx);
333:   /* Setup refinement */
334:   if (ctx->geometryCtx.areaCtx != PETSC_NULL) {
335:     MeshSetUserContext(mesh, geomCtx->areaCtx);
336:   } else {
337:     MeshSetUserContext(mesh, &geomCtx->maxArea);
338:   }
339:   /* Report on mesh */
340:   MeshGetPartition(mesh, &part);
341:   switch(ctx->dim) {
342:   case 1:
343:     PartitionGetTotalElements(part, &totalElements);
344:     PartitionGetTotalNodes(part, &totalNodes);
345:     PetscPrintf(comm, "Elements: %d Nodes: %d\n", totalElements, totalNodes);
346:     break;
347:   case 2:
348:     PartitionGetTotalElements(part, &totalElements);
349:     PartitionGetTotalNodes(part, &totalNodes);
350:     PartitionGetTotalEdges(part, &totalEdges);
351:     PetscPrintf(comm, "Elements: %d Nodes: %d Edges: %d\n", totalElements, totalNodes, totalEdges);
352:     break;
353:   default:
354:     SETERRQ1(PETSC_ERR_SUP, "Do not support meshes of dimension %d", ctx->dim);
355:   }

357:   *m = mesh;
358:   return(0);
359: }

361: #undef  __FUNCT__
363: /*@
364:   PoissonContextCreateGrid - This function creates a grid for the main problem.

366:   Collective on PoissonContext

368:   Input Parameter:
369: . ctx     - A PoissonContext with problem specific information

371:   Level: beginner

373: .seealso PoissonRefineGrid(), PoissonDestroyGrid()
374: @*/
375: int PoissonContextCreateGrid(PoissonContext ctx) {
376:   Mesh mesh;
377:   Grid grid;
378:   int  ierr;

381:   /* Construct the mesh */
382:   PoissonContextCreateMesh(ctx, &mesh);
383:   /* Construct the grid */
384:   GridCreate(mesh, &grid);
385:   MeshDestroy(mesh);
386:   switch(ctx->dim) {
387:   case 1:
388:     if (ctx->linear == PETSC_TRUE) {
389:       if (ctx->flipDisc == PETSC_FALSE) {
390:         GridAddField(grid, "u", DISCRETIZATION_TRIANGULAR_1D_CONSTANT, 1, PETSC_NULL);
391:         GridAddField(grid, "v", DISCRETIZATION_TRIANGULAR_1D_LINEAR,   1, PETSC_NULL);
392:       } else {
393:         GridAddField(grid, "u", DISCRETIZATION_TRIANGULAR_1D_LINEAR,   1, PETSC_NULL);
394:         GridAddField(grid, "v", DISCRETIZATION_TRIANGULAR_1D_CONSTANT, 1, PETSC_NULL);
395:       }
396:     } else {
397:       GridAddField(grid, "u", DISCRETIZATION_TRIANGULAR_1D_LINEAR,    1, PETSC_NULL);
398:       GridAddField(grid, "v", DISCRETIZATION_TRIANGULAR_1D_QUADRATIC, 1, PETSC_NULL);
399:     }
400:     break;
401:   case 2:
402:     if (ctx->linear == PETSC_TRUE) {
403:       GridAddField(grid, "u", DISCRETIZATION_TRIANGULAR_2D_CONSTANT, 2, PETSC_NULL);
404:       GridAddField(grid, "v", DISCRETIZATION_TRIANGULAR_2D_LINEAR,   2, PETSC_NULL);
405:     } else {
406:       GridAddField(grid, "u", DISCRETIZATION_TRIANGULAR_2D_LINEAR,    2, PETSC_NULL);
407:       GridAddField(grid, "v", DISCRETIZATION_TRIANGULAR_2D_QUADRATIC, 2, PETSC_NULL);
408:     }
409:     break;
410:   default:
411:     SETERRQ1(PETSC_ERR_SUP, "Do not support meshes of dimension %d", ctx->dim);
412:   }
413:   GridSetFromOptions(grid);

415:   ctx->grid = grid;
416:   return(0);
417: }

419: #undef  __FUNCT__
421: /*@
422:   PoissonRefineGrid - This function refines the mesh for the main grid.

424:   Collective on PoissonContext

426:   Input Parameters:
427: + ctx - A PoissonContext

429:   Options Database Keys:
430: . mesh_max_area - The maximum area na element may have

432:   Level: beginner

434: .seealso PoissonCreateGRid(), PoissonDestroyGrid()
435: @*/
436: int PoissonRefineGrid(PoissonContext ctx) {
437:   Grid                 oldGrid = ctx->grid;
438:   Grid                 grid;
439:   Mesh                 mesh;
440:   Partition            part;
441:   MeshGeometryContext *geomCtx = &ctx->geometryCtx;
442:   int                  totalElements, totalNodes, totalEdges;
443:   int                  ierr;

446:   /* Construct a refined mesh */
447:   if (geomCtx->areaCtx != PETSC_NULL) {
448:     GridRefineMesh(oldGrid, geomCtx->areaFunc, &grid);
449:     GridGetMesh(grid, &mesh);
450:     MeshSetUserContext(mesh, geomCtx->areaCtx);
451:   } else {
452:     geomCtx->maxArea *= 0.5;
453:     GridRefineMesh(oldGrid, geomCtx->areaFunc, &grid);
454:     GridGetMesh(grid, &mesh);
455:     MeshSetUserContext(mesh, &geomCtx->maxArea);
456:   }
457:   GridSetOptionsPrefix(grid, "ref_");
458:   GridSetFromOptions(grid);

460:   MeshGetPartition(mesh, &part);
461:   PartitionGetTotalElements(part, &totalElements);
462:   PartitionGetTotalNodes(part, &totalNodes);
463:   PartitionGetTotalEdges(part, &totalEdges);
464:   PetscPrintf(ctx->comm, "Elements: %d Nodes: %d Edges: %d\n", totalElements, totalNodes, totalEdges);
465:   CHKMEMQ;

467:   /* Replace old grid with refined grid */
468:   GridFinalizeBC(oldGrid);
469:   GridDestroy(oldGrid);
470:   ctx->grid = grid;
471:   return(0);
472: }

474: /*------------------------------------------------- Grid Setup -------------------------------------------------------*/
475: #undef  __FUNCT__
477: /*@
478:   PoissonSetupGrid - This function sets all the functions,
479:   operators , and boundary conditions for the problem.
480:   It also sets the parameters associated with the fields.

482:   Collective on Grid

484:   Input Parameters:
485: . grid - The grid
486: . ctx  - A PoissonContext

488:   Options Database Keys:
489: . use_laplacian - Use the Laplacian in stead of the Rate-of-Strain tensor

491:   Level: intermediate

493: .seealso PoissonCreateGrid()
494: @*/
495: int PoissonSetupGrid(PoissonContext ctx) {
496:   Grid grid   = ctx->grid;
497:   int  uField = 0;
498:   int  vField = 1;
499:   int  ierr;

502:   /* Setup Problem */
503:   GridSetActiveField(grid, uField);
504:   GridAddActiveField(grid, vField);

506:   /* Setup Rhs */
507:   PoissonSetupRhsFunction(grid, ctx);

509:   /* Setup Jacobian */
510:   if (ctx->flipDisc == PETSC_FALSE) {
511:     GridSetMatOperator(grid, vField, vField, IDENTITY,   -1.0, PETSC_FALSE, PETSC_NULL);
512:     GridAddMatOperator(grid, uField, vField, GRADIENT,   -1.0, PETSC_FALSE, PETSC_NULL);
513:     GridAddMatOperator(grid, vField, uField, DIVERGENCE,  1.0, PETSC_FALSE, PETSC_NULL);
514:   } else {
515:     GridSetMatOperator(grid, vField, vField, IDENTITY,   -1.0, PETSC_FALSE, PETSC_NULL);
516:     GridAddMatOperator(grid, uField, vField, DIVERGENCE, -1.0, PETSC_FALSE, PETSC_NULL);
517:     GridAddMatOperator(grid, vField, uField, GRADIENT,    1.0, PETSC_FALSE, PETSC_NULL);
518:   }

520:   /* Setup Dirchlet boundary conditions */
521:   PoissonSetupBC(grid, ctx);
522:   return(0);
523: }

525: #undef  __FUNCT__
527: /*@ PoissonSetupRhsFunction
528:   PoissonSetupRhsFunction - This function chooses a forcing  function for the problem.

530:   Collective on Grid

532:   Input Parameters:
533: . grid - The grid
534: . ctx  - A PoissonContext

536:   Level: intermediate

538:   Options Database Keys:

540: .seealso PoissonSetupGrid
541: @*/
542: int PoissonSetupRhsFunction(Grid grid, PoissonContext ctx) {
543:   int uField = 0;

547:   GridAddRhsFunction(grid, uField, RhsFunction, 1.0);
548:   return(0);
549: }

551: #undef  __FUNCT__
553: /*@ PoissonSetupBC
554:   PoissonSetupBC - This function chooses boundary conditions for the problem.

556:   Collective on Grid

558:   Input Parameters:
559: . grid - The grid
560: . ctx  - A PoissonContext

562:   Level: intermediate

564:   Options Database Keys:
565: . bc_reduce - Explicitly reduce the system using boundary conditions

567: .seealso PoissonSetupGrid()
568: @*/
569: int PoissonSetupBC(Grid grid, PoissonContext ctx) {
570:   int        uField = 0;
571:   int        vField = 1;
572:   PetscTruth reduceSystem;
573:   int        ierr;

576:   PetscOptionsHasName(PETSC_NULL, "-bc_reduce", &reduceSystem);
577:   switch(ctx->dim) {
578:   case 1:
579:     GridAddBC(grid, TOP_BD,    vField, GradientSolutionFunction, reduceSystem);
580:     GridAddBC(grid, BOTTOM_BD, vField, GradientSolutionFunction, reduceSystem);
581:     GridAddBC(grid, TOP_BD,    uField, SolutionFunction,         reduceSystem);
582:     GridAddBC(grid, BOTTOM_BD, uField, SolutionFunction,         reduceSystem);
583:     break;
584:   case 2:
585:     GridSetBC(grid, OUTER_BD, vField, GradientSolutionFunction, reduceSystem);
586:     break;
587:   default:
588:     SETERRQ1(PETSC_ERR_SUP, "Do not support meshes of dimension %d", ctx->dim);
589:   }
590:   GridSetBCContext(grid, ctx);
591:   return(0);
592: }

594: /*------------------------------------------------- SLES Setup -------------------------------------------------------*/
595: #undef  __FUNCT__
597: int PoissonCreateStructures(PoissonContext ctx) {

601:   /* Create the linear solver */
602:   SLESCreate(ctx->comm, &ctx->sles);
603:   SLESSetFromOptions(ctx->sles);
604:   /* Create solution, rhs, and exact solution vectors */
605:   GVecCreate(ctx->grid, &ctx->u);
606:   PetscObjectSetName((PetscObject) ctx->u, "Solution");
607:   GVecDuplicate(ctx->u, &ctx->f);
608:   PetscObjectSetName((PetscObject) ctx->f, "Rhs");
609:   GVecDuplicate(ctx->u, &ctx->uExact);
610:   PetscObjectSetName((PetscObject) ctx->uExact, "ExactSolution");
611:   /* Create the system matrix */
612:   GMatCreate(ctx->grid, &ctx->A);
613:   return(0);
614: }

616: #undef  __FUNCT__
618: int PoissonSetupKSP(KSP ksp, PoissonContext ctx) {
619:   GVecErrorKSPMonitorCtx *monCtx = &ctx->monitorCtx;
620:   PetscViewer             v;
621:   PetscDraw               draw;
622:   PetscTruth              opt;
623:   int                     ierr;

626:   /* Setup convergence monitors */
627:   PetscOptionsHasName(PETSC_NULL, "-error_viewer", &opt);
628:   if (opt == PETSC_TRUE) {
629:     v    = PETSC_VIEWER_DRAW_(ctx->comm);
630:     PetscViewerSetFormat(v, PETSC_VIEWER_DRAW_LG);
631:     PetscViewerDrawGetDraw(v, 0, &draw);
632:     PetscDrawSetTitle(draw, "Error");
633:   } else {
634:     v    = PETSC_VIEWER_STDOUT_(ctx->comm);
635:   }
636:   monCtx->error_viewer      = v;
637:   monCtx->solution          = ctx->uExact;
638:   monCtx->norm_error_viewer = PETSC_VIEWER_STDOUT_(ctx->comm);
639:   KSPSetMonitor(ksp, GVecErrorKSPMonitor, monCtx, PETSC_NULL);
640:   return(0);
641: }


644: #undef  __FUNCT__
646: int PoissonSetupPC(PC pc, PoissonContext ctx) {
647:   MPI_Comm    comm;
648:   PetscScalar alpha;
649:   int         uField = 0;
650:   int         ierr;

653:   PetscObjectGetComm((PetscObject) pc, &comm);
654:   GVecCreate(ctx->grid, &ctx->constantU);
655:   GVecEvaluateFunction(ctx->constantU, 1, &uField, PointFunctionOne, 1.0, PETSC_NULL);
656:   VecDot(ctx->constantU, ctx->constantU, &alpha);
657:   alpha = 1.0/PetscSqrtScalar(alpha);
658:   VecScale(&alpha, ctx->constantU);
659:   MatNullSpaceCreate(comm, PETSC_FALSE, 1, &ctx->constantU, &ctx->nullSpace);
660:   MatNullSpaceTest(ctx->nullSpace, ctx->A);
661:   PCNullSpaceAttach(pc, ctx->nullSpace);
662:   MatNullSpaceDestroy(ctx->nullSpace);
663:   return(0);
664: }

666: #undef  __FUNCT__
668: int PoissonSetupStructures(PoissonContext ctx) {
669:   KSP          ksp;
670:   PC           pc;
671:   int          uField = 0;
672:   int          vField = 1;
673:   MatStructure flag;
674:   PetscTruth   opt;
675:   int          ierr;

678:   /* Setup the linear solver */
679:   SLESSetOperators(ctx->sles, ctx->A, ctx->A, SAME_NONZERO_PATTERN);
680:   SLESGetKSP(ctx->sles, &ksp);
681:   PoissonSetupKSP(ksp, ctx);
682:   SLESGetPC(ctx->sles, &pc);
683:   PoissonSetupPC(pc, ctx);
684:   /* Evaluate the rhs */
685:   GridEvaluateRhs(ctx->grid, PETSC_NULL, ctx->f, (PetscObject) ctx);
686:   /* Evaluate the exact solution */
687:   GVecEvaluateFunction(ctx->uExact, 1, &uField, SolutionFunction,         1.0, ctx);
688:   GVecEvaluateFunction(ctx->uExact, 1, &vField, GradientSolutionFunction, 1.0, ctx);
689:   MatNullSpaceRemove(ctx->nullSpace, ctx->uExact, PETSC_NULL);
690:   /* Evaluate the system matrix */
691:   flag = DIFFERENT_NONZERO_PATTERN;
692:   GridEvaluateSystemMatrix(ctx->grid, PETSC_NULL, &ctx->A, &ctx->A, &flag, (PetscObject) ctx);
693:   MatCheckSymmetry(ctx->A);
694:   /* Apply Dirchlet boundary conditions */
695:   GMatSetBoundary(ctx->A, 1.0, ctx);
696:   GVecSetBoundary(ctx->f, ctx);
697:   /* View structures */
698:   PetscOptionsHasName(PETSC_NULL, "-mat_view", &opt);
699:   if (opt == PETSC_TRUE) {
700:     MatView(ctx->A, PETSC_VIEWER_STDOUT_(ctx->comm));
701:     MatView(ctx->A, PETSC_VIEWER_DRAW_(ctx->comm));
702:   }
703:   PetscOptionsHasName(PETSC_NULL, "-vec_view", &opt);
704:   if (opt == PETSC_TRUE) {
705:     VecView(ctx->f, PETSC_VIEWER_STDOUT_(ctx->comm));
706:   }
707:   return(0);
708: }

710: #undef  __FUNCT__
712: int PoissonDestroyStructures(PoissonContext ctx) {

716:   SLESDestroy(ctx->sles);
717:   GVecDestroy(ctx->f);
718:   GVecDestroy(ctx->u);
719:   GVecDestroy(ctx->uExact);
720:   GMatDestroy(ctx->A);
721:   GVecDestroy(ctx->constantU);
722:   return(0);
723: }

725: /*----------------------------------------------- Sanity Checks ------------------------------------------------------*/
726: #undef  __FUNCT__
728: int MatCheckSymmetry(Mat A) {
729:   Mat        trA;
730:   PetscTruth isSym;
731:   int        ierr;

734:   MatTranspose(A, &trA);
735:   MatEqual(A, trA, &isSym);
736:   MatDestroy(trA);
737:   if (isSym == PETSC_FALSE) PetscFunctionReturn(1);
738:   return(0);
739: }

741: #undef  __FUNCT__
743: int PoissonCheckSolution(PoissonContext ctx, GVec u, const char type[]) {
744:   GVec        r;
745:   PetscScalar minusOne = -1.0;
746:   PetscScalar norm;
747:   int         ierr;

750:   GVecDuplicate(u, &r);
751:   /* A u */
752:   MatMult(ctx->A, u, r);
753:   /* f - A u^* */
754:   VecAYPX(&minusOne, ctx->f, r);
755:   /* || f - A u || */
756:   VecNorm(r, NORM_2, &norm);
757:   PetscPrintf(ctx->comm, "Residual of the %s solution: %g\n", type, norm);

759:   GVecDestroy(r);
760:   return(0);
761: }

763: /*----------------------------------------------- Problem Callbacks --------------------------------------------------*/
764: #undef  __FUNCT__
766: /*@
767:   SolutionFunction - This function is the velocity solution function for the problem.

769:   Not collective

771:   Input Parameters:
772: + n      - The number of points
773: . comp   - The number of components
774: . x,y,z  - The points
775: . values - The output
776: - ctx    - A PoissonContext

778:   Level: beginner

780:   Note:
781: $  The solution for one dimension  is u = x^2 \hat x.
782: $  The solution for two dimensions is u = x^2 \hat x - 2 x y \hat y.

784: .keywords velocity, solution
785: .seealso RhsFunction()
786: @*/
787: int SolutionFunction(int n, int comp, double *x, double *y, double *z, PetscScalar *values, void *ctx) {
788:   PoissonContext pCtx = (PoissonContext) ctx;
789:   int            i;

792:   switch(comp) {
793:   case 1:
794:     if (pCtx->flipDisc == PETSC_FALSE) {
795:       for(i = 0; i < n; i++) {
796:         values[i] = x[i]*x[i];
797:       }
798:     } else {
799:       for(i = 0; i < n; i++) {
800:         values[i] = x[i]*x[i] - 1.5; // This is such a kludge, but I don't know how to do it automatically
801:       }
802:     }
803:     break;
804:   case 2:
805:     for(i = 0; i < n; i++) {
806:       values[i*2+0] = x[i]*x[i];
807:       values[i*2+1] = -2.0*x[i]*y[i];
808:     }
809:     break;
810:   default:
811:     SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid number of components %d", comp);
812:   }
813:   return(0);
814: }

816: #undef  __FUNCT__
818: /*@
819:   GradientSolutionFunction - This function is the solution function for the auxiliary field.

821:   Not collective

823:   Input Parameters:
824: + n      - The number of points
825: . comp   - The number of components
826: . x,y,z  - The points
827: . values - The output
828: - ctx    - A PoissonContext

830:   Level: beginner

832:   Note:
833: $  The solution for one dimension  is v = -2 x \hat x.
834: $  The solution for two dimensions is v = 

836: .keywords velocity, solution
837: .seealso RhsFunction()
838: @*/
839: int GradientSolutionFunction(int n, int comp, double *x, double *y, double *z, PetscScalar *values, void *ctx) {
840:   int i;

843:   switch(comp) {
844:   case 1:
845:     for(i = 0; i < n; i++) {
846:       values[i] = -2.0*x[i];
847:     }
848:     break;
849:   case 2:
850:     for(i = 0; i < n; i++) {
851:       values[i*2+0] = x[i]*x[i];
852:       values[i*2+1] = -2.0*x[i]*y[i];
853:     }
854:     break;
855:   default:
856:     SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid number of components %d", comp);
857:   }
858:   return(0);
859: }

861: #undef  __FUNCT__
863: /*@
864:   RhsFunction - This function is the forcing function for the problem.

866:   Not collective

868:   Input Parameters:
869: + n      - The number of points
870: . comp   - The number of components
871: . x,y,z  - The points
872: . values - The output
873: - ctx    - A PoissonContext

875:   Level: beginner

877:   Note:
878: $  The rhs for one dimension  is -\Delta u = -2 \hat x.
879: $  The rhs for two dimensions is -\Delta u = -2 \hat x.

881: .keywords velocity, rhs
882: .seealso SolutionFunction()
883: @*/
884: int RhsFunction(int n, int comp, double *x, double *y, double *z, PetscScalar *values, void *ctx) {
885:   /*PoissonContext s = (PoissonContext) ctx;*/
886:   int            i;

889:   switch(comp) {
890:   case 1:
891:     for(i = 0; i < n; i++) {
892:       values[i] = -2.0;
893:     }
894:     break;
895:   case 2:
896:     for(i = 0; i < n; i++) {
897:       values[i*2+0] = -2.0;
898:       values[i*2+1] = 0.0;
899:     }
900:     break;
901:   default:
902:     SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid number of components %d", comp);
903:   }
904:   return(0);
905: }

907: /*----------------------------------------------- Main Computation ---------------------------------------------------*/
908: #undef  __FUNCT__
910: int PoissonSolve(PoissonContext ctx, GVec f, GVec u, int *its) {

914:   /* Solve A u = f */
915:   SLESSolve(ctx->sles, f, u, its);

917:   /* Show solution */
918:   GVecViewFromOptions(ctx->uExact, "Exact Solution");
919:   GVecViewFromOptions(ctx->u,      "Solution");
920:   return(0);
921: }

923: #undef  __FUNCT__
925: int PoissonComputeField(PoissonContext ctx) {
926:   int its;  /* The iteration count for the linear solver */
927:   int loop;

930:   /* Get command-line options */
931:   PoissonContextSetup(ctx);

933:   for(loop = 0; loop < ctx->numLoops; loop++) {
934:     if (loop >= ctx->refineStep) {
935:       PoissonRefineGrid(ctx);
936:     }

938:     /* Setup problem */
939:     PoissonSetupGrid(ctx);
940:     PoissonCreateStructures(ctx);
941:     PoissonSetupStructures(ctx);

943:     /* Check the exact solution */
944:     PoissonCheckSolution(ctx, ctx->uExact, "exact");

946:     /* Solve system */
947:     PoissonSolve(ctx, ctx->f, ctx->u, &its);

949:     /* Check the computed solution */
950:     PoissonCheckSolution(ctx, ctx->u, "computed");

952:     /* Cleanup */
953:     PoissonDestroyStructures(ctx);
954:   }
955:   return(0);
956: }