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: }