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