Actual source code: poisson.c

  1: #ifdef PETSC_RCS_HEADER
  2: static char vcid[] = "$Id: poisson.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 poisson.h

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

 15:     -\Delta u = f(x)

 17:   Thus we now let

 19:     SolutionFunction() --> u_S (which we prescribe)
 20:     RhsFunction()      --> f
 21: */
 22: #undef  __FUNCT__
 24: int main(int argc, char **args) {
 25:   PoissonContext ctx; /* Holds problem specific information */
 26:   int           ierr;

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

 31:   PoissonInitializePackage(PETSC_NULL);                                      CHKERRABORT(PETSC_COMM_WORLD, ierr);
 32:   PoissonContextCreate(PETSC_COMM_WORLD, &ctx);                              CHKERRABORT(PETSC_COMM_WORLD, ierr);
 33:   PoissonComputeField(ctx);                                                  CHKERRABORT(PETSC_COMM_WORLD, ierr);
 34:   PoissonContextDestroy(ctx);                                                CHKERRABORT(PETSC_COMM_WORLD, ierr);

 36:   CHKMEMQ;
 37:   PetscFinalize();
 38:   return(0);
 39: }

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

 46:   Input Parameter:
 47: . path - The dynamic library path, or PETSC_NULL

 49:   Level: developer

 51: .keywords: Poisson, initialize, package
 52: .seealso: PetscInitialize()
 53: @*/
 54: int PoissonInitializePackage(char *path) {
 55:   static PetscTruth initialized = PETSC_FALSE;
 56:   char              logList[256];
 57:   char             *className;
 58:   PetscTruth        opt;
 59:   int               ierr;

 62:   if (initialized == PETSC_TRUE) return(0);
 63:   initialized = PETSC_TRUE;
 64:   /* Register Classes */
 65:   PetscLogClassRegister(&POISSON_COOKIE, "Poisson");
 66:   /* Register Constructors and Serializers */
 67:   /* Register Events */
 68:   PetscLogEventRegister(&POISSON_ComputeField, "PoissonComputeField", POISSON_COOKIE);
 69:   /* Process info exclusions */
 70:   PetscOptionsGetString(PETSC_NULL, "-log_info_exclude", logList, 256, &opt);
 71:   if (opt == PETSC_TRUE) {
 72:     PetscStrstr(logList, "PoissonContext", &className);
 73:     if (className) {
 74:       PetscLogInfoDeactivateClass(POISSON_COOKIE);
 75:     }
 76:   }
 77:   /* Process summary exclusions */
 78:   PetscOptionsGetString(PETSC_NULL, "-log_summary_exclude", logList, 256, &opt);
 79:   if (opt == PETSC_TRUE) {
 80:     PetscStrstr(logList, "PoissonContext", &className);
 81:     if (className) {
 82:       PetscLogEventDeactivateClass(POISSON_COOKIE);
 83:     }
 84:   }
 85:   return(0);
 86: }

 88: /*-------------------------------------------- PoissonContext Creation ------------------------------------------------*/
 89: #undef  __FUNCT__
 91: /*@
 92:   PoissonContextCreate - This function initializes the Poisson context.

 94:   Collective on MPI_Comm

 96:   Input Parameter:
 97: . comm - The communicator

 99:   Output Parameter:
100: . sCtx - The PoissonContext

102:   Level: beginner

104: .keywords: Poisson, context, create
105: .seealso: PoissonContextDestroy(), PoissonContextPrint(), PoissonContextSetup()
106: @*/
107: int PoissonContextCreate(MPI_Comm comm, PoissonContext *sCtx) {
108:   PoissonContext ctx;

111:   /* Setup context */
112:   PetscHeaderCreate(ctx, _PoissonContext, int, POISSON_COOKIE, -1, "PoissonContext", comm, 0, 0);
113:   PetscLogObjectCreate(ctx);
114:   PetscLogObjectMemory(ctx, sizeof(struct _PoissonContext));

116:   /* Initialize subobjects */
117:   ctx->grid       = PETSC_NULL;
118:   ctx->sles       = PETSC_NULL;
119:   ctx->A          = PETSC_NULL;
120:   ctx->u          = PETSC_NULL;
121:   ctx->f          = PETSC_NULL;
122:   ctx->uExact     = PETSC_NULL;
123:   /* Setup domain */
124:   ctx->geometryCtx.size[0]  = 2.0;
125:   ctx->geometryCtx.size[1]  = 2.0;
126:   ctx->geometryCtx.start[0] = 0.0;
127:   ctx->geometryCtx.start[1] = 0.0;
128:   /* Setup refinement */
129:   ctx->geometryCtx.maxArea  = 0.5;
130:   ctx->geometryCtx.areaFunc = PointFunctionConstant;
131:   ctx->geometryCtx.areaCtx  = PETSC_NULL;
132:   /* Initialize problem loop */
133:   ctx->dim        = 2;
134:   ctx->linear     = PETSC_FALSE;
135:   ctx->refineStep = 0;
136:   ctx->numLoops   = 0;

138:   *sCtx = ctx;
139:   return(0);
140: }

142: #undef  __FUNCT__
144: /*@
145:   PoissonContextDestroy - This function destroys the Poisson context.

147:   Collective on PoissonContext

149:   Input Parameter:
150: . ctx - The PoissonContext

152:   Level: beginner

154: .keywords: Poisson, context, destroy
155: .seealso: PoissonContextCreate(), PoissonContextSetup()
156: @*/
157: int PoissonContextDestroy(PoissonContext ctx) {
158:   Grid grid = ctx->grid;
159:   int  ierr;

162:   if (--ctx->refct > 0) SETERRQ(PETSC_ERR_PLIB, "Poisson context should not be referenced more than once");
163:   GridFinalizeBC(grid);
164:   GridDestroy(grid);
165:   PetscLogObjectDestroy(ctx);
166:   PetscHeaderDestroy(ctx);
167:   return(0);
168: }

170: /*--------------------------------------------- PoissonContext Setup -------------------------------------------------*/
171: #undef  __FUNCT__
173: /*@
174:   PoissonContextSetup - This function configures the context from user options.

176:   Collective on PoissonContext

178:   Input Parameter:
179: . ctx - The PoissonContext

181:   Options Database Keys:
182: + dim <num>           - The problem dimension
183: . linear              - Use a linear discretization
184: . num_systems <num>   - The number of systems to solve
185: . refine_step <num>   - The step to start refining the mesh
186: - mesh_max_area <num> - The maximum area of a triangle in the refined mesh

188:   Level: beginner

190: .seealso PoissonRefineGrid(), PoissonDestroyGrid()
191: */
192: int PoissonContextSetup(PoissonContext ctx) {
193:   PetscTruth opt;
194:   int        ierr;

197:   /* Determine the problem dimension */
198:   PetscOptionsGetInt(PETSC_NULL, "-dim", &ctx->dim, &opt);
199:   /* Determine the element type */
200:   PetscOptionsHasName(PETSC_NULL, "-linear", &ctx->linear);
201:   /* Determine how many systems to solve */
202:   PetscOptionsGetInt(PETSC_NULL, "-num_systems", &ctx->numLoops, &opt);
203:   /* The first iteration at which to refine the mesh */
204:   PetscOptionsGetInt(PETSC_NULL, "-refine_step", &ctx->refineStep, &opt);
205:   /* The maximum area of any triangle in the refined mesh */
206:   PetscOptionsGetReal("mesh_", "-max_area", &ctx->geometryCtx.maxArea, &opt);

208:   /* Create main problem */
209:   PoissonContextCreateGrid(ctx);
210:   return(0);
211: }

213: /*------------------------------------------------ Grid Creation -----------------------------------------------------*/
214: #undef  __FUNCT__
216: /*@
217:   PoissonContextCreateMeshBoundary - This function creates a mesh boundary for the main problem.

219:   Collective on PoissonContext

221:   Input Parameter:
222: . ctx - A PoissonContext with problem specific information

224:   Level: beginner

226: .seealso PoissonContextDestroyMeshBoundary(), PoissonContextCreateMesh()
227: @*/
228: int PoissonContextCreateMeshBoundary(PoissonContext ctx) {
229:   MPI_Comm             comm;
230:   MeshGeometryContext *geomCtx = &ctx->geometryCtx;
231:   int                  ierr;

234:   PetscObjectGetComm((PetscObject) ctx, &comm);
235:   switch(ctx->dim) {
236:   case 1:
237:     MeshBoundary1DCreateSimple(comm, geomCtx, &ctx->boundaryCtx);
238:     break;
239:   case 2:
240:     MeshBoundary2DCreateSimple(comm, geomCtx, &ctx->boundaryCtx);
241:     break;
242:   default:
243:     SETERRQ1(PETSC_ERR_SUP, "Do not support meshes of dimension %d", ctx->dim);
244:   }
245:   return(0);
246: }

248: #undef  __FUNCT__
250: /*@
251:   PoissonContextDestroyMeshBoundary - This function destroys the mesh boundary for the main problem.

253:   Collective on PoissonContext

255:   Input Parameter:
256: . ctx - A PoissonContext with problem specific information

258:   Level: beginner

260: .seealso PoissonContextCreateMeshBoundary(), PoissonContextCreateMesh()
261: @*/
262: int PoissonContextDestroyMeshBoundary(PoissonContext ctx) {
263:   MPI_Comm comm;
264:   int      ierr;

267:   PetscObjectGetComm((PetscObject) ctx, &comm);
268:   switch(ctx->dim) {
269:   case 1:
270:     MeshBoundary1DDestroy(comm, &ctx->boundaryCtx);
271:     break;
272:   case 2:
273:     MeshBoundary2DDestroy(comm, &ctx->boundaryCtx);
274:     break;
275:   default:
276:     SETERRQ1(PETSC_ERR_SUP, "Do not support meshes of dimension %d", ctx->dim);
277:   }
278:   return(0);
279: }

281: #undef  __FUNCT__
283: /*@
284:   PoissonContextCreateMesh - This function creates a mesh for the main problem.

286:   Collective on PoissonContext

288:   Input Parameter:
289: . ctx - A PoissonContext with problem specific information

291:   Output Parameter:
292: . m   - The Mesh

294:   Options Database Keys:
295: . mesh_refine   - Refines the mesh based on area criteria
296: . mesh_max_area - The maximum area of an element

298:   Level: beginner

300: .seealso PoissonRefineGrid(), PoissonDestroyGrid()
301: @*/
302: int PoissonContextCreateMesh(PoissonContext ctx, Mesh *m) {
303:   MeshGeometryContext *geomCtx = &ctx->geometryCtx;
304:   MPI_Comm             comm;
305:   Mesh                 mesh;
306:   Partition            part;
307:   int                  totalElements, totalNodes, totalEdges;
308:   char                 name[1024];
309:   int                  d;
310:   int                  ierr;

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

367:   *m = mesh;
368:   return(0);
369: }

371: #undef  __FUNCT__
373: /*@
374:   PoissonContextCreateGrid - This function creates a grid for the main problem.

376:   Collective on PoissonContext

378:   Input Parameter:
379: . ctx     - A PoissonContext with problem specific information

381:   Level: beginner

383: .seealso PoissonRefineGrid(), PoissonDestroyGrid()
384: @*/
385: int PoissonContextCreateGrid(PoissonContext ctx) {
386:   Mesh mesh;
387:   Grid grid;
388:   int  ierr;

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

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

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

425:   Collective on PoissonContext

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

430:   Options Database Keys:
431: . mesh_max_area - The maximum area an element may have

433:   Level: beginner

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

448:   /* Construct a refined mesh */
449:   if (geomCtx->areaCtx != PETSC_NULL) {
450:     GridRefineMesh(oldGrid, geomCtx->areaFunc, &grid);
451:     GridGetMesh(grid, &mesh);
452:     MeshSetUserContext(mesh, geomCtx->areaCtx);
453:   } else {
454:     geomCtx->maxArea *= 0.5;
455:     GridRefineMesh(oldGrid, geomCtx->areaFunc, &grid);
456:     GridGetMesh(grid, &mesh);
457:     MeshSetUserContext(mesh, &geomCtx->maxArea);
458:   }
459:   sprintf(name, "mesh.r%.4g", geomCtx->maxArea);
460:   PetscObjectSetName((PetscObject) mesh, name);
461:   MeshSetOptionsPrefix(mesh, "ref_");
462:   GridSetOptionsPrefix(grid, "ref_");
463:   GridSetFromOptions(grid);

465:   MeshGetPartition(mesh, &part);
466:   PartitionGetTotalElements(part, &totalElements);
467:   PartitionGetTotalNodes(part, &totalNodes);
468:   PartitionGetTotalEdges(part, &totalEdges);
469:   PetscPrintf(ctx->comm, "Elements: %d Nodes: %d Edges: %d\n", totalElements, totalNodes, totalEdges);
470:   CHKMEMQ;

472:   /* Replace old grid with refined grid */
473:   GridFinalizeBC(oldGrid);
474:   GridDestroy(oldGrid);
475:   ctx->grid = grid;
476:   return(0);
477: }

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

487:   Collective on Grid

489:   Input Parameters:
490: + grid - The grid
491: - ctx  - A PoissonContext

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

496:   Level: intermediate

498: .seealso PoissonCreateGrid()
499: @*/
500: int PoissonSetupGrid(PoissonContext ctx) {
501:   Grid grid  = ctx->grid;
502:   int  field = 0;
503:   int  ierr;

506:   /* Setup Problem */
507:   GridSetActiveField(grid, field);

509:   /* Setup Rhs */
510:   PoissonSetupRhsFunction(grid, ctx);

512:   /* Setup Jacobian */
513:   GridSetMatOperator(grid, field, field, LAPLACIAN, -1.0, PETSC_FALSE, PETSC_NULL);

515:   /* Setup Dirchlet boundary conditions */
516:   PoissonSetupBC(grid, ctx);
517:   return(0);
518: }

520: #undef  __FUNCT__
522: /*@ PoissonSetupRhsFunction
523:   PoissonSetupRhsFunction - This function chooses a forcing  function for the problem.

525:   Collective on Grid

527:   Input Parameters:
528: + grid - The grid
529: - ctx  - A PoissonContext

531:   Level: intermediate

533:   Options Database Keys:

535: .seealso PoissonSetupGrid
536: @*/
537: int PoissonSetupRhsFunction(Grid grid, PoissonContext ctx) {
538:   int field = 0;

542:   GridAddRhsFunction(grid, field, RhsFunction, 1.0);
543:   return(0);
544: }

546: #undef  __FUNCT__
548: /*@ PoissonSetupBC
549:   PoissonSetupBC - This function chooses boundary conditions for the problem.

551:   Collective on Grid

553:   Input Parameters:
554: + grid - The grid
555: - ctx  - A PoissonContext

557:   Level: intermediate

559:   Options Database Keys:
560: . bc_reduce - Explicitly reduce the system using boundary conditions

562: .seealso PoissonSetupGrid()
563: @*/
564: int PoissonSetupBC(Grid grid, PoissonContext ctx) {
565:   int        field   = 0;
566:   PetscTruth reduceSystem;
567:   int        ierr;

570:   PetscOptionsHasName(PETSC_NULL, "-bc_reduce", &reduceSystem);
571:   switch(ctx->dim) {
572:   case 1:
573:     GridSetBC(grid, TOP_BD,    field, SolutionFunction, reduceSystem);
574:     GridAddBC(grid, BOTTOM_BD, field, SolutionFunction, reduceSystem);
575:     break;
576:   case 2:
577:     GridSetBC(grid, OUTER_BD, field, SolutionFunction, reduceSystem);
578:     break;
579:   default:
580:     SETERRQ1(PETSC_ERR_SUP, "Do not support meshes of dimension %d", ctx->dim);
581:   }
582:   GridSetBCContext(grid, ctx);
583:   return(0);
584: }

586: /*------------------------------------------------- SLES Setup -------------------------------------------------------*/
587: #undef  __FUNCT__
589: int PoissonCreateStructures(PoissonContext ctx) {

593:   /* Create the linear solver */
594:   SLESCreate(ctx->comm, &ctx->sles);
595:   SLESSetFromOptions(ctx->sles);
596:   /* Create solution, rhs, and exact solution vectors */
597:   GVecCreate(ctx->grid, &ctx->u);
598:   PetscObjectSetName((PetscObject) ctx->u, "Solution");
599:   GVecDuplicate(ctx->u, &ctx->f);
600:   PetscObjectSetName((PetscObject) ctx->f, "Rhs");
601:   GVecDuplicate(ctx->u, &ctx->uExact);
602:   PetscObjectSetName((PetscObject) ctx->uExact, "ExactSolution");
603:   /* Create the system matrix */
604:   GMatCreate(ctx->grid, &ctx->A);
605:   return(0);
606: }

608: #undef  __FUNCT__
610: int PoissonSetupKSP(KSP ksp, PoissonContext ctx) {
611:   GVecErrorKSPMonitorCtx *monCtx = &ctx->monitorCtx;
612:   PetscViewer             v;
613:   PetscDraw               draw;
614:   PetscTruth              opt;
615:   int                     ierr;

618:   /* Setup convergence monitors */
619:   PetscOptionsHasName(PETSC_NULL, "-error_viewer", &opt);
620:   if (opt == PETSC_TRUE) {
621:     v    = PETSC_VIEWER_DRAW_(ctx->comm);
622:     PetscViewerSetFormat(v, PETSC_VIEWER_DRAW_LG);
623:     PetscViewerDrawGetDraw(v, 0, &draw);
624:     PetscDrawSetTitle(draw, "Error");
625:   } else {
626:     v    = PETSC_VIEWER_STDOUT_(ctx->comm);
627:   }
628:   monCtx->error_viewer      = v;
629:   monCtx->solution          = ctx->uExact;
630:   monCtx->norm_error_viewer = PETSC_VIEWER_STDOUT_(ctx->comm);
631:   KSPSetMonitor(ksp, GVecErrorKSPMonitor, monCtx, PETSC_NULL);
632:   return(0);
633: }


636: #undef  __FUNCT__
638: int PoissonSetupPC(PC pc, PoissonContext ctx) {
640:   return(0);
641: }

643: #undef  __FUNCT__
645: int PoissonSetupStructures(PoissonContext ctx) {
646:   KSP          ksp;
647:   PC           pc;
648:   int          field  = 0;
649:   MatStructure flag;
650:   PetscTruth   opt;
651:   int          ierr;

654:   /* Setup the linear solver */
655:   SLESSetOperators(ctx->sles, ctx->A, ctx->A, SAME_NONZERO_PATTERN);
656:   SLESGetKSP(ctx->sles, &ksp);
657:   PoissonSetupKSP(ksp, ctx);
658:   SLESGetPC(ctx->sles, &pc);
659:   PoissonSetupPC(pc, ctx);
660:   /* Evaluate the rhs */
661:   GridEvaluateRhs(ctx->grid, PETSC_NULL, ctx->f, (PetscObject) ctx);
662:   /* Evaluate the exact solution */
663:   GVecEvaluateFunction(ctx->uExact, 1, &field, SolutionFunction, 1.0, ctx);
664:   /* Evaluate the system matrix */
665:   flag = DIFFERENT_NONZERO_PATTERN;
666:   GridEvaluateSystemMatrix(ctx->grid, PETSC_NULL, &ctx->A, &ctx->A, &flag, (PetscObject) ctx);
667:   MatCheckSymmetry(ctx->A);
668:   /* Apply Dirchlet boundary conditions */
669:   GMatSetBoundary(ctx->A, 1.0, ctx);
670:   GVecSetBoundary(ctx->f, ctx);
671:   /* View structures */
672:   PetscOptionsHasName(PETSC_NULL, "-mat_view", &opt);
673:   if (opt == PETSC_TRUE) {
674:     MatView(ctx->A, PETSC_VIEWER_STDOUT_(ctx->comm));
675:     MatView(ctx->A, PETSC_VIEWER_DRAW_(ctx->comm));
676:   }
677:   PetscOptionsHasName(PETSC_NULL, "-vec_view", &opt);
678:   if (opt == PETSC_TRUE) {
679:     VecView(ctx->f, PETSC_VIEWER_STDOUT_(ctx->comm));
680:   }
681:   PetscOptionsHasName(PETSC_NULL, "-system_view", &opt);
682:   if (opt == PETSC_TRUE) {
683:     PetscViewer viewer;
684:     char        filename[256];
685:     int         len;

687:     PetscOptionsGetString(PETSC_NULL, "-system_view", filename, 255, &opt);
688:     PetscStrlen(filename, &len);
689:     if (len > 0) {
690:       PetscViewerBinaryOpen(ctx->comm, filename, PETSC_BINARY_CREATE, &viewer);
691:     } else {
692:       PetscViewerBinaryOpen(ctx->comm, "system.petsc", PETSC_BINARY_CREATE, &viewer);
693:     }
694:     MatView(ctx->A, viewer);
695:     VecView(ctx->f, viewer);
696:     PetscViewerDestroy(viewer);
697:   }
698:   return(0);
699: }

701: #undef  __FUNCT__
703: int PoissonDestroyStructures(PoissonContext ctx) {

707:   SLESDestroy(ctx->sles);
708:   GVecDestroy(ctx->f);
709:   GVecDestroy(ctx->u);
710:   GVecDestroy(ctx->uExact);
711:   GMatDestroy(ctx->A);
712:   return(0);
713: }

715: /*----------------------------------------------- Sanity Checks ------------------------------------------------------*/
716: #undef  __FUNCT__
718: int MatCheckSymmetry(Mat A) {
719:   Mat        trA;
720:   PetscTruth isSym;
721:   int        ierr;

724:   MatTranspose(A, &trA);
725:   MatEqual(A, trA, &isSym);
726:   MatDestroy(trA);
727:   if (isSym == PETSC_FALSE) PetscFunctionReturn(1);
728:   return(0);
729: }

731: #undef  __FUNCT__
733: int PoissonCheckSolution(PoissonContext ctx, GVec u, const char type[]) {
734:   GVec        r;
735:   PetscScalar minusOne = -1.0;
736:   PetscScalar norm;
737:   int         ierr;

740:   GVecDuplicate(u, &r);
741:   /* A u */
742:   MatMult(ctx->A, u, r);
743:   /* f - A u^* */
744:   VecAYPX(&minusOne, ctx->f, r);
745:   /* || f - A u || */
746:   VecNorm(r, NORM_2, &norm);
747:   PetscPrintf(ctx->comm, "Residual of the %s solution: %g\n", type, norm);

749:   GVecDestroy(r);
750:   return(0);
751: }

753: /*----------------------------------------------- Problem Callbacks --------------------------------------------------*/

755: #undef  __FUNCT__
757: /*@
758:   SolutionFunction - This function is the velocity solution function for the problem.

760:   Not collective

762:   Input Parameters:
763: + n      - The number of points
764: . comp   - The number of components
765: . x,y,z  - The points
766: . values - The output
767: - ctx    - A PoissonContext

769:   Level: beginner

771:   Note:
772: $  The solution for one dimension  is u = x^2 \hat x.
773: $  The solution for two dimensions is u = x^2 \hat x - 2 x y \hat y.

775: .keywords velocity, solution
776: .seealso RhsFunction()
777: @*/
778: int SolutionFunction(int n, int comp, double *x, double *y, double *z, PetscScalar *values, void *ctx) {
779:   int i;

782:   switch(comp) {
783:   case 1:
784:     for(i = 0; i < n; i++) {
785:       values[i] = x[i]*x[i];
786:     }
787:     break;
788:   case 2:
789:     for(i = 0; i < n; i++) {
790:       values[i*2+0] = x[i]*x[i];
791:       values[i*2+1] = -2.0*x[i]*y[i];
792:     }
793:     break;
794:   default:
795:     SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid number of components %d", comp);
796:   }
797:   return(0);
798: }

800: #undef  __FUNCT__
802: /*@
803:   RhsFunction - This function is the forcing function for the problem.

805:   Not collective

807:   Input Parameters:
808: + n      - The number of points
809: . comp   - The number of components
810: . x,y,z  - The points
811: . values - The output
812: - ctx    - A PoissonContext

814:   Level: beginner

816:   Note:
817: $  The rhs for one dimension  is -\Delta u = -2 \hat x.
818: $  The rhs for two dimensions is -\Delta u = -2 \hat x.

820: .keywords velocity, rhs
821: .seealso SolutionFunction()
822: @*/
823: int RhsFunction(int n, int comp, double *x, double *y, double *z, PetscScalar *values, void *ctx) {
824:   /*PoissonContext s = (PoissonContext) ctx;*/
825:   int            i;

828:   switch(comp) {
829:   case 1:
830:     for(i = 0; i < n; i++) {
831:       values[i] = -2.0;
832:     }
833:     break;
834:   case 2:
835:     for(i = 0; i < n; i++) {
836:       values[i*2+0] = -2.0;
837:       values[i*2+1] = 0.0;
838:     }
839:     break;
840:   default:
841:     SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid number of components %d", comp);
842:   }
843:   return(0);
844: }

846: /*----------------------------------------------- Main Computation ---------------------------------------------------*/
847: #undef  __FUNCT__
849: int PoissonSolve(PoissonContext ctx, GVec f, GVec u, int *its) {

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

856:   /* Show solution */
857:   GVecViewFromOptions(ctx->uExact, "Exact Solution");
858:   GVecViewFromOptions(ctx->u,      "Solution");
859:   return(0);
860: }

862: #undef  __FUNCT__
864: int PoissonComputeField(PoissonContext ctx) {
865:   int its;  /* The iteration count for the linear solver */
866:   int loop;

869:   /* Get command-line options */
870:   PoissonContextSetup(ctx);

872:   for(loop = 0; loop < ctx->numLoops; loop++) {
873:     if (loop >= ctx->refineStep) {
874:       PoissonRefineGrid(ctx);
875:     }

877:     /* Setup problem */
878:     PoissonSetupGrid(ctx);
879:     PoissonCreateStructures(ctx);
880:     PoissonSetupStructures(ctx);

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

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

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

891:     /* Cleanup */
892:     PoissonDestroyStructures(ctx);
893:   }
894:   return(0);
895: }