Actual source code: stokes.c
1: #ifdef PETSC_RCS_HEADER
2: static char vcid[] = "$Id: stokes.c,v 1.5 2000/02/01 17:02:51 knepley Exp $";
3: #endif
5: static char help[] = "This is the Stokes problem with P2/P1 elements.\n\n";
7: int STOKES_COOKIE;
8: int STOKES_ComputeField;
10: #include stokes.h
12: /*
13: Here we are solving the Stokes equation:
15: / -\nu \Delta u - \nabla p \ = / f(x)\
16: \ \nabla \cdot u / \ 0 /
18: Thus we now let
20: VelocitySolutionFunction() --> u_S (which we prescribe)
21: PressureSolutionFunction() --> p_S (which we prescribe)
22: VelocityRhsFunction() --> f
23: */
24: #undef __FUNCT__
26: int main(int argc, char **args) {
27: StokesContext ctx; /* Holds problem specific information */
28: int ierr;
31: PetscInitialize(&argc, &args, 0, help); CHKERRABORT(PETSC_COMM_WORLD, ierr);
33: StokesInitializePackage(PETSC_NULL); CHKERRABORT(PETSC_COMM_WORLD, ierr);
34: StokesContextCreate(PETSC_COMM_WORLD, &ctx); CHKERRABORT(PETSC_COMM_WORLD, ierr);
35: StokesComputeFlowField(ctx); CHKERRABORT(PETSC_COMM_WORLD, ierr);
36: StokesContextDestroy(ctx); CHKERRABORT(PETSC_COMM_WORLD, ierr);
38: CHKMEMQ;
39: PetscFinalize();
40: return(0);
41: }
45: /*@C
46: StokesInitializePackage - This function initializes everything in the Stokes package.
48: Input Parameter:
49: . path - The dynamic library path, or PETSC_NULL
51: Level: developer
53: .keywords: Stokes, initialize, package
54: .seealso: PetscInitialize()
55: @*/
56: int StokesInitializePackage(char *path) {
57: static PetscTruth initialized = PETSC_FALSE;
58: char logList[256];
59: char *className;
60: PetscTruth opt;
61: int ierr;
64: if (initialized == PETSC_TRUE) return(0);
65: initialized = PETSC_TRUE;
66: /* Register Classes */
67: PetscLogClassRegister(&STOKES_COOKIE, "Stokes");
68: /* Register Constructors and Serializers */
69: /* Register Events */
70: PetscLogEventRegister(&STOKES_ComputeField, "StokesComputeField", STOKES_COOKIE);
71: /* Process info exclusions */
72: PetscOptionsGetString(PETSC_NULL, "-log_info_exclude", logList, 256, &opt);
73: if (opt == PETSC_TRUE) {
74: PetscStrstr(logList, "StokesContext", &className);
75: if (className) {
76: PetscLogInfoDeactivateClass(STOKES_COOKIE);
77: }
78: }
79: /* Process summary exclusions */
80: PetscOptionsGetString(PETSC_NULL, "-log_summary_exclude", logList, 256, &opt);
81: if (opt == PETSC_TRUE) {
82: PetscStrstr(logList, "StokesContext", &className);
83: if (className) {
84: PetscLogEventDeactivateClass(STOKES_COOKIE);
85: }
86: }
87: return(0);
88: }
90: /*-------------------------------------------- StokesContext Creation ------------------------------------------------*/
91: #undef __FUNCT__
93: /*@
94: StokesContextCreate - This function initializes the Stokes context.
96: Collective on MPI_Comm
98: Input Parameter:
99: . comm - The communicator
101: Output Parameter:
102: . sCtx - The StokesContext
104: Level: beginner
106: .keywords: Stokes, context, create
107: .seealso: StokesContextDestroy(), StokesContextPrint(), StokesContextSetup()
108: @*/
109: int StokesContextCreate(MPI_Comm comm, StokesContext *sCtx) {
110: StokesContext ctx;
113: /* Setup context */
114: PetscHeaderCreate(ctx, _StokesContext, int, STOKES_COOKIE, -1, "context", comm, 0, 0);
115: PetscLogObjectCreate(ctx);
116: PetscLogObjectMemory(ctx, sizeof(struct _StokesContext));
118: /* Initialize subobjects */
119: ctx->grid = PETSC_NULL;
120: ctx->sles = PETSC_NULL;
121: ctx->A = PETSC_NULL;
122: ctx->u = PETSC_NULL;
123: ctx->p = PETSC_NULL;
124: ctx->f = PETSC_NULL;
125: ctx->origF = PETSC_NULL;
126: ctx->uExact = PETSC_NULL;
127: ctx->pExact = PETSC_NULL;
128: ctx->RofS = -1;
129: /* Setup domain */
130: ctx->geometryCtx.size[0] = 2.0;
131: ctx->geometryCtx.size[1] = 2.0;
132: ctx->geometryCtx.start[0] = 0.0;
133: ctx->geometryCtx.start[1] = 0.0;
134: /* Setup refinement */
135: ctx->geometryCtx.maxArea = 0.5;
136: ctx->geometryCtx.areaFunc = PointFunctionConstant;
137: ctx->geometryCtx.areaCtx = PETSC_NULL;
138: /* Initialize physical information */
139: ctx->physicalCtx.nu = -1.0;
140: ctx->physicalCtx.fluidSG = 1.0;
141: /* Initialize preconditioning suites */
142: ctx->useML = PETSC_FALSE;
143: ctx->useSchur = PETSC_FALSE;
144: /* Initialize problem loop */
145: ctx->dim = 2;
146: ctx->linear = PETSC_FALSE;
147: ctx->refineStep = 0;
148: ctx->numLoops = 0;
150: *sCtx = ctx;
151: return(0);
152: }
154: #undef __FUNCT__
156: /*@
157: StokesContextDestroy - This function destroys the Stokes context.
159: Collective on StokesContext
161: Input Parameter:
162: . ctx - The StokesContext
164: Level: beginner
166: .keywords: Stokes, context, destroy
167: .seealso: StokesContextCreate(), StokesContextSetup()
168: @*/
169: int StokesContextDestroy(StokesContext ctx) {
170: Grid grid = ctx->grid;
171: int ierr;
174: if (--ctx->refct > 0) SETERRQ(PETSC_ERR_PLIB, "Stokes context should not be referenced more than once");
175: GridFinalizeBC(grid);
176: GridDestroy(grid);
177: PetscLogObjectDestroy(ctx);
178: PetscHeaderDestroy(ctx);
179: return(0);
180: }
182: /*--------------------------------------------- StokesContext Setup --------------------------------------------------*/
183: #undef __FUNCT__
185: /*@
186: StokesContextSetup - This function configures the context from user options.
188: Collective on StokesContext
190: Input Parameter:
191: . ctx - A StokesContext
193: Options Database Keys:
194: + dim <num> - The problem dimension
195: . linear - Using a linear discretization
196: . num_systems <num> - The number of systems to solve
197: . refine_step <num> - The step to start refining the mesh
198: . physical_nu <num> - The fluid kinematic viscosity
199: . physical_fluid_sg <num> - The fluid specific gravity
200: . mesh_max_area <num> - The maximum area of a triangle in the refined mesh
201: . pc_ml - The ML preconditioning suite
202: - pc_schur - The Schur preconditioning suite
204: Level: beginner
206: .seealso StokesRefineGrid(), StokesDestroyGrid()
207: @*/
208: int StokesContextSetup(StokesContext ctx) {
209: PetscTruth opt;
210: int ierr;
213: /* Determine the problem dimension */
214: PetscOptionsGetInt(PETSC_NULL, "-dim", &ctx->dim, &opt);
215: /* Determine the element type */
216: PetscOptionsHasName(PETSC_NULL, "-linear", &ctx->linear);
217: /* The first iteration at which to refine the mesh */
218: PetscOptionsGetInt(PETSC_NULL, "-refine_step", &ctx->refineStep, &opt);
219: /* Determine how many systems to solve */
220: PetscOptionsGetInt(PETSC_NULL, "-num_systems", &ctx->numLoops, &opt);
221: /* Setup refinement */
222: PetscOptionsGetReal("mesh", "-max_area", &ctx->geometryCtx.maxArea, &opt);
223: /* Get physical information */
224: PetscOptionsGetReal("physical_", "-nu", &ctx->physicalCtx.nu, &opt);
225: PetscOptionsGetReal("physical_", "-fluid_sg", &ctx->physicalCtx.fluidSG, &opt);
226: /* Setup preconditioning suites */
227: PetscOptionsHasName(PETSC_NULL, "-pc_ml", &ctx->useML);
228: PetscOptionsHasName(PETSC_NULL, "-pc_schur", &ctx->useSchur);
230: /* Create main problem */
231: StokesContextCreateGrid(ctx);
232: return(0);
233: }
235: /*------------------------------------------------ Grid Creation -----------------------------------------------------*/
236: #undef __FUNCT__
238: /*@
239: StokesContextCreateMeshBoundary - This function creates a mesh boundary for the main problem.
241: Collective on StokesContext
243: Input Parameter:
244: . ctx - The StokesContext with problem specific information
246: Level: beginner
248: .seealso StokesContextDestroyMeshBoundary(), StokesContextCreateMesh()
249: @*/
250: int StokesContextCreateMeshBoundary(StokesContext ctx) {
251: MPI_Comm comm;
252: MeshGeometryContext *geomCtx = &ctx->geometryCtx;
253: int ierr;
256: PetscObjectGetComm((PetscObject) ctx, &comm);
257: switch(ctx->dim) {
258: case 1:
259: MeshBoundary1DCreateSimple(comm, geomCtx, &ctx->boundaryCtx);
260: break;
261: case 2:
262: MeshBoundary2DCreateSimple(comm, geomCtx, &ctx->boundaryCtx);
263: break;
264: default:
265: SETERRQ1(PETSC_ERR_SUP, "Do not support meshes of dimension %d", ctx->dim);
266: }
267: return(0);
268: }
270: #undef __FUNCT__
272: /*@
273: StokesContextDestroyMeshBoundary - This function destroys a mesh boundary for the main problem.
275: Collective on StokesContext
277: Input Parameter:
278: . ctx - The StokesContext with problem specific information
280: Level: beginner
282: .seealso StokesContextCreateMeshBoundary(), StokesContextDestroyMesh()
283: @*/
284: int StokesContextDestroyMeshBoundary(StokesContext ctx) {
285: MPI_Comm comm;
286: int ierr;
289: PetscObjectGetComm((PetscObject) ctx, &comm);
290: switch(ctx->dim) {
291: case 1:
292: MeshBoundary1DDestroy(comm, &ctx->boundaryCtx);
293: break;
294: case 2:
295: MeshBoundary2DDestroy(comm, &ctx->boundaryCtx);
296: break;
297: default:
298: SETERRQ1(PETSC_ERR_SUP, "Do not support meshes of dimension %d", ctx->dim);
299: }
300: return(0);
301: }
303: #undef __FUNCT__
305: /*@
306: StokesContextCreateMesh - This function creates a mesh for the main problem.
308: Collective on StokesContext
310: Input Parameter:
311: . ctx - A StokesContext with problem specific information
313: Output Parameter:
314: . m - The Mesh
316: Options Database Keys:
317: . mesh_refine - Refines the mesh based on area criteria
318: . mesh_max_area - The maximum area of an element
320: Level: beginner
322: .seealso StokesRefineGrid(), StokesDestroyGrid()
323: @*/
324: int StokesContextCreateMesh(StokesContext ctx, Mesh *m) {
325: MeshGeometryContext *geomCtx = &ctx->geometryCtx;
326: MPI_Comm comm;
327: Mesh mesh;
328: Partition part;
329: int totalElements, totalNodes, totalEdges;
330: char name[1024];
331: int d;
332: int ierr;
335: PetscObjectGetComm((PetscObject) ctx, &comm);
336: StokesContextCreateMeshBoundary(ctx);
337: MeshCreate(comm, &mesh);
338: MeshSetDimension(mesh, ctx->dim);
339: for(d = 0; d < ctx->dim; d++) {
340: MeshSetPeriodicDimension(mesh, d, geomCtx->isPeriodic[d]);
341: }
342: switch(ctx->dim) {
343: case 1:
344: if (ctx->linear == PETSC_TRUE) {
345: MeshSetNumCorners(mesh, 2);
346: } else {
347: MeshSetNumCorners(mesh, 3);
348: }
349: break;
350: case 2:
351: if (ctx->linear == PETSC_TRUE) {
352: MeshSetNumCorners(mesh, 3);
353: } else {
354: MeshSetNumCorners(mesh, 6);
355: }
356: break;
357: default:
358: SETERRQ1(PETSC_ERR_SUP, "Do not support meshes of dimension %d", ctx->dim);
359: }
360: MeshSetBoundary(mesh, &ctx->boundaryCtx);
361: sprintf(name, "mesh.r%.4g", geomCtx->maxArea);
362: PetscObjectSetName((PetscObject) mesh, name);
363: MeshSetFromOptions(mesh);
364: StokesContextDestroyMeshBoundary(ctx);
365: /* Setup refinement */
366: if (ctx->geometryCtx.areaCtx != PETSC_NULL) {
367: MeshSetUserContext(mesh, geomCtx->areaCtx);
368: } else {
369: MeshSetUserContext(mesh, &geomCtx->maxArea);
370: }
371: /* Report on mesh */
372: MeshGetPartition(mesh, &part);
373: switch(ctx->dim) {
374: case 1:
375: PartitionGetTotalElements(part, &totalElements);
376: PartitionGetTotalNodes(part, &totalNodes);
377: PetscPrintf(ctx->comm, "Elements: %d Nodes: %d\n", totalElements, totalNodes);
378: break;
379: case 2:
380: PartitionGetTotalElements(part, &totalElements);
381: PartitionGetTotalNodes(part, &totalNodes);
382: PartitionGetTotalEdges(part, &totalEdges);
383: PetscPrintf(ctx->comm, "Elements: %d Nodes: %d Edges: %d\n", totalElements, totalNodes, totalEdges);
384: break;
385: default:
386: SETERRQ1(PETSC_ERR_SUP, "Do not support meshes of dimension %d", ctx->dim);
387: }
389: *m = mesh;
390: return(0);
391: }
393: #undef __FUNCT__
395: /*@
396: StokesContextCreateGrid - This function creates a grid for the main problem.
398: Collective on StokesContext
400: Input Parameter:
401: . ctx - A StokesContext with problem specific information
403: Level: beginner
405: .seealso StokesRefineGrid(), StokesDestroyGrid()
406: @*/
407: int StokesContextCreateGrid(StokesContext ctx) {
408: Mesh mesh;
409: Grid grid;
410: int ierr;
413: /* Construct the mesh */
414: StokesContextCreateMesh(ctx, &mesh);
415: /* Construct the grid */
416: GridCreate(mesh, &grid);
417: MeshDestroy(mesh);
419: switch(ctx->dim) {
420: case 1:
421: if (ctx->linear == PETSC_TRUE) {
422: GridAddField(grid, "Velocity", DISCRETIZATION_TRIANGULAR_1D_LINEAR, 2, PETSC_NULL);
423: GridAddField(grid, "Pressure", DISCRETIZATION_TRIANGULAR_1D_CONSTANT, 1, PETSC_NULL);
424: } else {
425: GridAddField(grid, "Velocity", DISCRETIZATION_TRIANGULAR_1D_QUADRATIC, 2, PETSC_NULL);
426: GridAddField(grid, "Pressure", DISCRETIZATION_TRIANGULAR_1D_LINEAR, 1, PETSC_NULL);
427: }
428: break;
429: case 2:
430: if (ctx->linear == PETSC_TRUE) {
431: GridAddField(grid, "Velocity", DISCRETIZATION_TRIANGULAR_2D_LINEAR, 2, PETSC_NULL);
432: GridAddField(grid, "Pressure", DISCRETIZATION_TRIANGULAR_2D_CONSTANT, 1, PETSC_NULL);
433: } else {
434: GridAddField(grid, "Velocity", DISCRETIZATION_TRIANGULAR_2D_QUADRATIC, 2, PETSC_NULL);
435: GridAddField(grid, "Pressure", DISCRETIZATION_TRIANGULAR_2D_LINEAR, 1, PETSC_NULL);
436: }
437: break;
438: default:
439: SETERRQ1(PETSC_ERR_SUP, "Do not support meshes of dimension %d", ctx->dim);
440: }
441: GridSetFromOptions(grid);
443: ctx->grid = grid;
444: return(0);
445: }
447: #undef __FUNCT__
449: /*@
450: StokesRefineGrid - This function refines the mesh for the main grid.
452: Collective on StokesContext
454: Input Parameters:
455: + ctx - A StokesContext
457: Options Database Keys:
458: . mesh_max_area - The maximum area an element may have
460: Level: beginner
462: .seealso StokesCreateGrid(), StokesDestroyGrid()
463: @*/
464: int StokesRefineGrid(StokesContext ctx) {
465: Grid oldGrid = ctx->grid;
466: Grid grid;
467: Mesh mesh;
468: Partition part;
469: MeshGeometryContext *geomCtx = &ctx->geometryCtx;
470: char name[1024];
471: int totalElements, totalNodes, totalEdges;
472: int ierr;
475: /* Construct a refined mesh */
476: if (geomCtx->areaCtx != PETSC_NULL) {
477: GridRefineMesh(oldGrid, geomCtx->areaFunc, &grid);
478: GridGetMesh(grid, &mesh);
479: MeshSetUserContext(mesh, geomCtx->areaCtx);
480: } else {
481: geomCtx->maxArea *= 0.5;
482: GridRefineMesh(oldGrid, geomCtx->areaFunc, &grid);
483: GridGetMesh(grid, &mesh);
484: MeshSetUserContext(mesh, &geomCtx->maxArea);
485: }
486: sprintf(name, "mesh.r%.4g", geomCtx->maxArea);
487: PetscObjectSetName((PetscObject) mesh, name);
488: MeshSetOptionsPrefix(mesh, "ref_");
489: GridSetOptionsPrefix(grid, "ref_");
490: GridSetFromOptions(grid);
492: MeshGetPartition(mesh, &part);
493: PartitionGetTotalElements(part, &totalElements);
494: PartitionGetTotalNodes(part, &totalNodes);
495: PartitionGetTotalEdges(part, &totalEdges);
496: PetscPrintf(ctx->comm, "Elements: %d Nodes: %d Edges: %d\n", totalElements, totalNodes, totalEdges);
497: CHKMEMQ;
499: /* Replace old grid with refined grid */
500: GridFinalizeBC(oldGrid);
501: GridDestroy(oldGrid);
502: ctx->grid = grid;
503: return(0);
504: }
506: /*------------------------------------------------- Grid Setup -------------------------------------------------------*/
507: #undef __FUNCT__
509: /*@
510: StokesSetupGrid - This function sets all the functions, operators , and boundary conditions for the problem.
511: It also sets the parameters associated with the fields.
513: Collective on Grid
515: Input Parameters:
516: + grid - The grid
517: - ctx - A StokesContext
519: Options Database Keys:
520: . use_laplacian - Use the Laplacian instead of the Rate-of-Strain tensor
522: Level: intermediate
524: .seealso StokesCreateGrid()
525: @*/
526: int StokesSetupGrid(StokesContext ctx) {
527: Grid grid = ctx->grid;
528: int vel = 0;
529: int pres = 1;
530: double nu = ctx->physicalCtx.nu;
531: double invRho = 1.0/ctx->physicalCtx.fluidSG;
532: PetscTruth opt;
533: int ierr;
536: /* Setup Fields */
537: GridSetFieldName(grid, vel, "Velocity");
538: GridSetFieldName(grid, pres, "Pressure");
540: /* Setup Problem */
541: GridSetActiveField(grid, vel);
542: if (ctx->useML == PETSC_FALSE) {
543: GridAddActiveField(grid, pres);
544: }
546: /* Setup Rhs */
547: PetscOptionsHasName(PETSC_NULL, "-use_laplacian", &opt);
548: if (opt) {
549: ctx->RofS = LAPLACIAN;
550: } else {
551: GridRegisterOperator(grid, vel, RateOfStrainTensor, &ctx->RofS);
552: }
554: GridSetRhsOperator(grid, vel, vel, ctx->RofS, -nu, PETSC_FALSE, PETSC_NULL);
555: if (ctx->useML == PETSC_FALSE) {
556: GridAddRhsOperator(grid, pres, vel, GRADIENT, -invRho, PETSC_FALSE, PETSC_NULL);
557: GridAddRhsOperator(grid, vel, pres, DIVERGENCE, invRho, PETSC_FALSE, PETSC_NULL);
558: }
559: StokesSetupRhsFunction(grid, ctx);
561: /* Setup Jacobian */
562: GridSetMatOperator(grid, vel, vel, ctx->RofS, -nu, PETSC_FALSE, PETSC_NULL);
563: if (ctx->useML == PETSC_FALSE) {
564: GridAddMatOperator(grid, pres, vel, GRADIENT, -invRho, PETSC_FALSE, PETSC_NULL);
565: GridAddMatOperator(grid, vel, pres, DIVERGENCE, invRho, PETSC_FALSE, PETSC_NULL);
566: }
568: /* Setup Dirchlet boundary conditions */
569: StokesSetupBC(grid, ctx);
570: return(0);
571: }
573: #undef __FUNCT__
575: /*@ StokesSetupRhsFunction
576: StokesSetupRhsFunction - This function chooses a forcing function for the problem.
578: Collective on Grid
580: Input Parameters:
581: + grid - The grid
582: - ctx - A StokesContext
584: Level: intermediate
586: Options Database Keys:
588: .seealso StokesSetupGrid
589: @*/
590: int StokesSetupRhsFunction(Grid grid, StokesContext ctx) {
591: int vel = 0;
592: int pres = 1;
596: GridAddRhsFunction(grid, vel, VelocityRhsFunction, 1.0);
597: if (ctx->useML == PETSC_FALSE) {
598: GridAddRhsFunction(grid, pres, PointFunctionZero, 1.0);
599: }
600: return(0);
601: }
603: #undef __FUNCT__
605: /*@ StokesSetupBC
606: StokesSetupBC - This function chooses boundary conditions for the problem.
608: Collective on Grid
610: Input Parameters:
611: + grid - The grid
612: - ctx - A StokesContext
614: Level: intermediate
616: Options Database Keys:
617: + bc_reduce - Explicitly reduce the system using boundary conditions
618: . bc_exact - Use the exact solution for boundary conditions
619: . bc_pressure_none - Do not impose boundary conditions on the pressure
620: . bc_pressure_full - Impose boundary conditions over the entire boundary
621: . bc_pressure_exact - Use the exact solution for boundary conditions
622: . bc_pressure_x - The x-coordinate of the pressure normalization point
623: . bc_pressure_y - The y-coordinate of the pressure normalization point
624: - bc_pressure_z - The z-coordinate of the pressure normalization point
626: .seealso StokesSetupGrid()
627: @*/
628: int StokesSetupBC(Grid grid, StokesContext ctx) {
629: MeshGeometryContext *geomCtx = &ctx->geometryCtx;
630: int vel = 0;
631: int pres = 1;
632: PetscTruth reduceSystem;
633: double pX, pY, pZ;
634: PetscTruth opt;
635: int ierr;
638: /* Setup reduction usng boundary conditions */
639: PetscOptionsHasName(PETSC_NULL, "-bc_reduce", &reduceSystem);
641: PetscOptionsHasName(PETSC_NULL, "-bc_exact", &opt);
642: if (opt == PETSC_TRUE) {
643: GridSetBC(grid, OUTER_BD, vel, VelocitySolutionFunction, reduceSystem);
644: } else {
645: GridSetBC(grid, OUTER_BD, vel, PointFunctionZero, reduceSystem);
646: }
648: /* Pressure boundary conditions */
649: PetscOptionsHasName(PETSC_NULL, "-bc_pressure_none", &opt);
650: if (opt == PETSC_FALSE) {
651: PetscOptionsHasName(PETSC_NULL, "-bc_pressure_full", &opt);
652: if (opt == PETSC_TRUE) {
653: GridAddBC(grid, OUTER_BD, pres, PressureSolutionFunction, reduceSystem);
654: } else {
655: pX = geomCtx->start[0];
656: pY = geomCtx->start[1];
657: pZ = geomCtx->start[2];
658: PetscOptionsGetReal(PETSC_NULL, "-bc_pressure_x", &pX, &opt);
659: PetscOptionsGetReal(PETSC_NULL, "-bc_pressure_y", &pY, &opt);
660: PetscOptionsGetReal(PETSC_NULL, "-bc_pressure_z", &pZ, &opt);
661: PetscOptionsHasName(PETSC_NULL, "-bc_pressure_exact", &opt);
662: if (opt == PETSC_TRUE) {
663: GridSetPointBC(grid, pX, pY, pZ, pres, PressureSolutionFunction, reduceSystem);
664: } else {
665: GridSetPointBC(grid, pX, pY, pZ, pres, PointFunctionZero, reduceSystem);
666: }
667: }
668: }
669: GridSetBCContext(grid, ctx);
670: return(0);
671: }
673: /*------------------------------------------------- SLES Setup -------------------------------------------------------*/
674: #undef __FUNCT__
676: int StokesCreateStructures(StokesContext ctx) {
677: VarOrdering pOrder;
678: int pres = 1;
679: int ierr;
682: /* Create the linear solver */
683: SLESCreate(ctx->comm, &ctx->sles);
684: SLESSetFromOptions(ctx->sles);
685: /* Create solution, rhs, projected rhs, and exact solution vectors */
686: GVecCreate(ctx->grid, &ctx->u);
687: PetscObjectSetName((PetscObject) ctx->u, "Solution");
688: GVecDuplicate(ctx->u, &ctx->f);
689: PetscObjectSetName((PetscObject) ctx->f, "Rhs");
690: GVecDuplicate(ctx->u, &ctx->origF);
691: PetscObjectSetName((PetscObject) ctx->origF, "OriginalRhs");
692: GVecDuplicate(ctx->u, &ctx->uExact);
693: PetscObjectSetName((PetscObject) ctx->uExact, "ExactSolution");
694: if (ctx->useML == PETSC_TRUE) {
695: FieldClassMap cm, reduceCM;
696: VarOrdering reduceOrder;
697: Vec reduceVec;
699: /* Create the multiplier, and the exact multiplier */
700: VarOrderingCreateGeneral(ctx->grid, 1, &pres, &pOrder);
701: GVecCreateRectangular(ctx->grid, pOrder, &ctx->p);
702: PetscObjectSetName((PetscObject) ctx->p, "Pressure");
703: GVecDuplicate(ctx->p, &ctx->pExact);
704: PetscObjectSetName((PetscObject) ctx->pExact, "ExactPressure");
706: VarOrderingGetClassMap(pOrder, &cm);
707: FieldClassMapReduce(cm, ctx->grid, &reduceCM);
708: VarOrderingCreateReduceGeneral(ctx->grid, cm, reduceCM, &reduceOrder);
709: GVecCreateRectangularGhost(ctx->grid, reduceOrder, &reduceVec);
710: GridCalcBCValues_Private(ctx->grid, reduceOrder, reduceVec, PETSC_FALSE, PETSC_NULL);
711: PetscObjectCompose((PetscObject) ctx->p, "reductionOrder", (PetscObject) reduceOrder);
712: PetscObjectCompose((PetscObject) ctx->pExact, "reductionOrder", (PetscObject) reduceOrder);
713: PetscObjectCompose((PetscObject) ctx->p, "reductionVec", (PetscObject) reduceVec);
714: PetscObjectCompose((PetscObject) ctx->pExact, "reductionVec", (PetscObject) reduceVec);
715: FieldClassMapDestroy(reduceCM);
716: VarOrderingDestroy(pOrder);
717: VarOrderingDestroy(reduceOrder);
718: VecDestroy(reduceVec);
719: }
720: /* Create the system matrix */
721: GMatCreate(ctx->grid, &ctx->A);
722: return(0);
723: }
725: #undef __FUNCT__
727: int StokesSetupKSP(KSP ksp, StokesContext ctx) {
728: GVecErrorKSPMonitorCtx *monCtx = &ctx->monitorCtx;
729: PetscViewer v;
730: PetscDraw draw;
731: PetscTruth opt;
732: int ierr;
735: /* Setup the multilevel preconditioner */
736: if (ctx->useML == PETSC_TRUE) {
737: KSPSetPreconditionerSide(ksp, PC_SYMMETRIC);
738: }
739: /* Setup convergence monitors */
740: PetscOptionsHasName(PETSC_NULL, "-error_viewer", &opt);
741: if (opt == PETSC_TRUE) {
742: v = PETSC_VIEWER_DRAW_(ctx->comm);
743: PetscViewerSetFormat(v, PETSC_VIEWER_DRAW_LG);
744: PetscViewerDrawGetDraw(v, 0, &draw);
745: PetscDrawSetTitle(draw, "Error");
746: } else {
747: v = PETSC_VIEWER_STDOUT_(ctx->comm);
748: }
749: monCtx->error_viewer = v;
750: monCtx->solution = ctx->uExact;
751: monCtx->norm_error_viewer = PETSC_VIEWER_STDOUT_(ctx->comm);
752: KSPSetMonitor(ksp, GVecErrorKSPMonitor, monCtx, PETSC_NULL);
753: return(0);
754: }
757: #undef __FUNCT__
759: int StokesSetupPC(PC pc, StokesContext ctx) {
760: int vel = 0;
761: int pres = 1;
762: PetscTruth opt;
763: int ierr;
766: /* Setup the multilevel preconditioner */
767: if (ctx->useML == PETSC_TRUE) {
768: #ifndef PETSC_USE_DYANMIC_LIBRARIES
769: GSolverInitializePackage(PETSC_NULL);
770: #endif
771: PetscOptionsHasName(PETSC_NULL, "-bc_reduce", &opt);
772: if (opt == PETSC_FALSE) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "The -bc_reduce option is necessary for ML");
774: PCSetType(pc, PCMULTILEVEL);
775: PCMultiLevelSetFields(pc, pres, vel);
776: PCMultiLevelSetGradientOperator(pc, GRADIENT, DIVERGENCE, -1.0);
777: PCSetVector(pc, ctx->f);
778: }
779: if (ctx->useSchur == PETSC_TRUE) {
780: PCSetType(pc, PCSCHUR);
781: PCSchurSetGradientOperator(pc, GRADIENT, DIVERGENCE);
782: }
783: return(0);
784: }
786: #undef __FUNCT__
788: int StokesSetupStructures(StokesContext ctx) {
789: KSP ksp;
790: PC pc;
791: VarOrdering pOrder;
792: int vel = 0;
793: int pres = 1;
794: MatStructure flag;
795: PetscTruth opt;
796: int ierr;
799: /* Setup the linear solver */
800: SLESSetOperators(ctx->sles, ctx->A, ctx->A, SAME_NONZERO_PATTERN);
801: SLESGetKSP(ctx->sles, &ksp);
802: StokesSetupKSP(ksp, ctx);
803: SLESGetPC(ctx->sles, &pc);
804: StokesSetupPC(pc, ctx);
805: /* Evaluate the rhs */
806: GridEvaluateRhs(ctx->grid, PETSC_NULL, ctx->f, (PetscObject) ctx);
807: VecCopy(ctx->f, ctx->origF);
808: /* Evaluate the exact solution */
809: GVecEvaluateFunction(ctx->uExact, 1, &vel, VelocitySolutionFunction, 1.0, ctx);
810: if (ctx->useML == PETSC_FALSE) {
811: GVecEvaluateFunction(ctx->uExact, 1, &pres, PressureSolutionFunction, 1.0, ctx);
812: }
813: /* Evaluate the exact multiplier */
814: if (ctx->useML == PETSC_TRUE) {
815: VarOrderingCreateGeneral(ctx->grid, 1, &pres, &pOrder);
816: GVecEvaluateFunctionRectangular(ctx->pExact, pOrder, PressureSolutionFunction, 1.0, ctx);
817: VarOrderingDestroy(pOrder);
818: }
819: /* Evaluate the system matrix */
820: flag = DIFFERENT_NONZERO_PATTERN;
821: GridEvaluateSystemMatrix(ctx->grid, PETSC_NULL, &ctx->A, &ctx->A, &flag, (PetscObject) ctx);
822: MatCheckSymmetry(ctx->A);
823: /* Apply Dirchlet boundary conditions */
824: GMatSetBoundary(ctx->A, 1.0, ctx);
825: GVecSetBoundary(ctx->f, ctx);
826: /* View structures */
827: PetscOptionsHasName(PETSC_NULL, "-mat_view", &opt);
828: if (opt == PETSC_TRUE) {
829: MatView(ctx->A, PETSC_VIEWER_STDOUT_(ctx->comm));
830: MatView(ctx->A, PETSC_VIEWER_DRAW_(ctx->comm));
831: }
832: PetscOptionsHasName(PETSC_NULL, "-vec_view", &opt);
833: if (opt == PETSC_TRUE) {
834: VecView(ctx->f, PETSC_VIEWER_STDOUT_(ctx->comm));
835: }
836: return(0);
837: }
839: #undef __FUNCT__
841: int StokesDestroyStructures(StokesContext ctx) {
845: SLESDestroy(ctx->sles);
846: GVecDestroy(ctx->f);
847: GVecDestroy(ctx->origF);
848: GVecDestroy(ctx->u);
849: GVecDestroy(ctx->uExact);
850: if (ctx->useML == PETSC_TRUE) {
851: GVecDestroy(ctx->p);
852: GVecDestroy(ctx->pExact);
853: }
854: GMatDestroy(ctx->A);
855: return(0);
856: }
858: /*----------------------------------------------- Sanity Checks ------------------------------------------------------*/
859: #undef __FUNCT__
861: int MatCheckSymmetry(Mat A) {
862: Mat trA;
863: PetscTruth isSym;
864: int ierr;
867: MatTranspose(A, &trA);
868: MatEqual(A, trA, &isSym);
869: MatDestroy(trA);
870: if (isSym == PETSC_FALSE) PetscFunctionReturn(1);
871: return(0);
872: }
874: #undef __FUNCT__
876: int StokesCheckSolution(StokesContext ctx, GVec u, GVec p, const char type[]) {
877: GVec r, projR;
878: PC pc;
879: PetscScalar minusOne = -1.0;
880: PetscScalar norm;
881: int ierr;
884: GVecDuplicate(u, &r);
885: GVecDuplicate(u, &projR);
886: SLESGetPC(ctx->sles, &pc);
888: /* Check the velocity solution */
889: /* A u^* */
890: MatMult(ctx->A, u, r);
891: /* f - A u^* */
892: VecAYPX(&minusOne, ctx->origF, r);
893: /* P^T_2 (f - A u^*) */
894: if (ctx->useML == PETSC_TRUE) {
895: PCApplySymmetricLeft(pc, r, projR);
896: VecNorm(projR, NORM_2, &norm);
897: PetscPrintf(ctx->comm, "Residual of the %s velocity solution: %lf\n", type, norm);
898: } else {
899: VecNorm(r, NORM_2, &norm);
900: PetscPrintf(ctx->comm, "Residual of the %s solution: %lf\n", type, norm);
901: }
903: /* Check the multiplier */
904: if (ctx->useML == PETSC_TRUE) {
905: /* P^T_2 B p^* */
906: PCMultiLevelApplyGradient(pc, p, r);
907: PCApplySymmetricLeft(pc, r, projR);
908: VecNorm(projR, NORM_2, &norm);
909: PetscPrintf(ctx->comm, "Residual of the %s gradient contribution: %lf\n", type, norm);
910: /* f - A u^* - B p^* */
911: MatMult(ctx->A, u, projR);
912: VecAYPX(&minusOne, ctx->origF, projR);
913: VecAXPY(&minusOne, r, projR);
914: VecNorm(projR, NORM_2, &norm);
915: PetscPrintf(ctx->comm, "Residual of the full %s solution: %lf\n", type, norm);
916: }
917:
918: GVecDestroy(r);
919: GVecDestroy(projR);
920: return(0);
921: }
923: /*----------------------------------------------- Problem Callbacks --------------------------------------------------*/
925: #undef __FUNCT__
927: /*@
928: VelocitySolutionFunction - This function is the velocity solution function for the problem.
930: Not collective
932: Input Parameters:
933: + n - The number of points
934: . comp - The number of components
935: . x,y,z - The points
936: . values - The output
937: - ctx - A StokesContext
939: Level: beginner
941: Note:
942: The solution is u = x^2 \hat x - 2 x y \hat y
944: .keywords velocity, solution
945: .seealso PressureSolutionFunction(), VelocityRhsFunction()
946: @*/
947: int VelocitySolutionFunction(int n, int comp, double *x, double *y, double *z, PetscScalar *values, void *ctx) {
948: int i;
951: if (comp != 2) SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid number of components %d", comp);
952: for(i = 0; i < n; i++) {
953: values[i*2+0] = x[i]*x[i];
954: values[i*2+1] = -2.0*x[i]*y[i];
955: }
956: return(0);
957: }
959: #undef __FUNCT__
961: /*@
962: PressureSolutionFunction - This function is the pressure solution function for the problem.
964: Not collective
966: Input Parameters:
967: + n - The number of points
968: . comp - The number of components
969: . x,y,z - The points
970: . values - The output
971: - ctx - A StokesContext
973: Level: beginner
975: Note:
976: The solution is p = x + y - 2, the integral of p is zero over the domain.
978: .keywords pressure, solution
979: .seealso VelocitySolutionFunction(), VelocityRhsFunction()
980: @*/
981: int PressureSolutionFunction(int n, int comp, double *x, double *y, double *z, PetscScalar *values, void *ctx) {
982: int i;
985: if (comp != 1) SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid number of components %d", comp);
986: for(i = 0; i < n; i++) {
987: values[i] = x[i] + y[i] - 2.0;
988: }
989: return(0);
990: }
992: #undef __FUNCT__
994: /*@
995: VelocityRhsFunction - This function is the forcing function for the problem.
997: Not collective
999: Input Parameters:
1000: + n - The number of points
1001: . comp - The number of components
1002: . x,y,z - The points
1003: . values - The output
1004: - ctx - A StokesContext
1006: Level: beginner
1008: Note:
1009: The rhs is -\nu \Delta u - \rho^{-1} \nabla p = -(2\nu + 1) \hat x - \hat y
1011: .keywords velocity, rhs
1012: .seealso VelocitySolutionFunction(), PressureSolutionFunction()
1013: @*/
1014: int VelocityRhsFunction(int n, int comp, double *x, double *y, double *z, PetscScalar *values, void *ctx) {
1015: StokesContext s = (StokesContext) ctx;
1016: double nu = s->physicalCtx.nu;
1017: double invRho = 1.0/s->physicalCtx.fluidSG;
1018: int i;
1021: if (comp != 2) SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid number of components %d", comp);
1022: for (i = 0; i < n; i++) {
1023: values[i*2+0] = -2.0*nu - invRho;
1024: values[i*2+1] = -invRho;
1025: }
1026: return(0);
1027: }
1029: /*-------------------------------------------------- Operators -------------------------------------------------------*/
1030: /* This defines the linear operator
1032: {1 \over 2} \nabla \cdot ( \nabla u + {\nabla u}^T )
1033: */
1034: #undef __FUNCT__
1036: int RateOfStrainTensor(Discretization disc, Discretization test, int rowSize, int colSize,
1037: int globalRowStart, int globalColStart, int globalSize, double *coords,
1038: PetscScalar alpha, PetscScalar *field, PetscScalar *array, void *ctx)
1039: {
1040: int comp; /* The number of field components */
1041: int funcs; /* The number of shape functions */
1042: int numQuadPoints; /* Number of points used for Gaussian quadrature */
1043: double *quadWeights; /* Weights in the standard element for Gaussian quadrature */
1044: double *quadShapeFuncDers; /* Shape function derivatives evaluated at quadrature points */
1045: double *quadTestFuncDers; /* Test function derivatives evaluated at quadrature points */
1046: double dxxi; /* \PartDer{x}{\xi} */
1047: double dxet; /* \PartDer{x}{\eta} */
1048: double dyxi; /* \PartDer{y}{\xi} */
1049: double dyet; /* \PartDer{y}{\eta} */
1050: double dxix; /* \PartDer{\xi}{x} */
1051: double detx; /* \PartDer{\eta}{x} */
1052: double dxiy; /* \PartDer{\xi}{y} */
1053: double dety; /* \PartDer{\eta}{y} */
1054: double dphix; /* \PartDer{\phi_j}{x} Shape */
1055: double dphiy; /* \PartDer{\phi_j}{y} Shape */
1056: double dpsix; /* \PartDer{\psi_i}{x} Test */
1057: double dpsiy; /* \PartDer{\psi_i}{y} Test */
1058: double jac; /* |J| for map to standard element */
1059: double invjac; /* |J^{-1}| for map from standard element */
1060: int i, j, f, p;
1061: int ierr;
1064: DiscretizationGetNumComponents(disc, &comp);
1065: DiscretizationGetNumFunctions(disc, &funcs);
1066: if (comp != 2) SETERRQ(PETSC_ERR_ARG_WRONG, "Operator only valid with 2D vector field");
1068: DiscretizationGetNumQuadraturePoints(disc, &numQuadPoints);
1069: DiscretizationGetQuadratureWeights(disc, &quadWeights);
1070: DiscretizationGetQuadratureDerivatives(disc, &quadShapeFuncDers);
1071: DiscretizationGetQuadratureDerivatives(test, &quadTestFuncDers);
1073: /* Bring out a factor of 1/2 */
1074: alpha *= 0.5;
1076: /* Calculate element matrix entries by Gaussian quadrature */
1077: for(p = 0; p < numQuadPoints; p++) {
1078: /* \PartDer{x}{\xi}(p) = \sum^{funcs}_{f=1} x_f \PartDer{\phi^f(p)}{\xi}
1079: \PartDer{x}{\eta}(p) = \sum^{funcs}_{f=1} x_f \PartDer{\phi^f(p)}{\eta}
1080: \PartDer{y}{\xi}(p) = \sum^{funcs}_{f=1} y_f \PartDer{\phi^f(p)}{\xi}
1081: \PartDer{y}{\eta}(p) = \sum^{funcs}_{f=1} y_f \PartDer{\phi^f(p)}{\eta} */
1082: dxxi = 0.0; dxet = 0.0;
1083: dyxi = 0.0; dyet = 0.0;
1084: for(f = 0; f < funcs; f++) {
1085: dxxi += coords[f*2] *quadShapeFuncDers[p*funcs*2+f*2];
1086: dxet += coords[f*2] *quadShapeFuncDers[p*funcs*2+f*2+1];
1087: dyxi += coords[f*2+1]*quadShapeFuncDers[p*funcs*2+f*2];
1088: dyet += coords[f*2+1]*quadShapeFuncDers[p*funcs*2+f*2+1];
1089: }
1090: jac = fabs(dxxi*dyet - dxet*dyxi);
1091: #ifdef PETSC_USE_BOPT_g
1092: if (jac < 1.0e-14) {
1093: PetscPrintf(PETSC_COMM_SELF, "p: %d x1: %lf y1: %lf x2: %lf y2: %lf x3: %lf y3: %lf\n",
1094: p, coords[0], coords[1], coords[2], coords[3], coords[4], coords[5]);
1095: SETERRQ(PETSC_ERR_DISC_SING_JAC, "Singular Jacobian");
1096: }
1097: #endif
1098: /* These are the elements of the inverse matrix */
1099: invjac = 1.0/jac;
1100: dxix = dyet*invjac;
1101: dxiy = -dxet*invjac;
1102: detx = -dyxi*invjac;
1103: dety = dxxi*invjac;
1105: for(i = 0; i < funcs; i++) {
1106: for(j = 0; j < funcs; j++) {
1107: dpsix = quadTestFuncDers[p*funcs*2+i*2] *dxix + quadTestFuncDers[p*funcs*2+i*2+1] *detx;
1108: dphix = quadShapeFuncDers[p*funcs*2+j*2]*dxix + quadShapeFuncDers[p*funcs*2+j*2+1]*detx;
1109: dpsiy = quadTestFuncDers[p*funcs*2+i*2] *dxiy + quadTestFuncDers[p*funcs*2+i*2+1] *dety;
1110: dphiy = quadShapeFuncDers[p*funcs*2+j*2]*dxiy + quadShapeFuncDers[p*funcs*2+j*2+1]*dety;
1111: array[(i*comp+0+globalRowStart)*globalSize + j*comp+0+globalColStart] +=
1112: -alpha*(2.0*dpsix*dphix + dpsiy*dphiy)*jac*quadWeights[p];
1113: array[(i*comp+0+globalRowStart)*globalSize + j*comp+1+globalColStart] +=
1114: -alpha*( dpsiy*dphix )*jac*quadWeights[p];
1115: array[(i*comp+1+globalRowStart)*globalSize + j*comp+0+globalColStart] +=
1116: -alpha*( dpsix*dphiy)*jac*quadWeights[p];
1117: array[(i*comp+1+globalRowStart)*globalSize + j*comp+1+globalColStart] +=
1118: -alpha*( dpsix*dphix + 2.0*dpsiy*dphiy)*jac*quadWeights[p];
1119: }
1120: }
1121: }
1122: PetscLogFlops((8*funcs + 8 + 38*funcs*funcs) * numQuadPoints);
1124: return(0);
1125: }
1127: /*----------------------------------------------- Main Computation ---------------------------------------------------*/
1128: #undef __FUNCT__
1130: int StokesSolve(StokesContext ctx, GVec f, GVec u, int *its) {
1131: PC pc;
1135: /* Solve P^T_2 A P_2 (P^{-1} u)_2 = P^T_2 (f - A P_1 D^{-1} Z^T g) */
1136: SLESSolve(ctx->sles, f, u, its);
1138: if (ctx->useML == PETSC_TRUE) {
1139: /* Calculate p = Z D^{-1} P^T_1 (f - A P_2 (P^{-1} u)_2 - A P_1 D^{-1} Z^T g) */
1140: SLESGetPC(ctx->sles, &pc);
1141: PCMultiLevelGetMultiplier(pc, u, ctx->p);
1143: /* Recover u */
1144: PCMultiLevelBuildSolution(pc, u);
1146: /* Show pressure */
1147: GVecViewFromOptions(ctx->pExact, "Exact Pressure");
1148: GVecViewFromOptions(ctx->p, "Pressure");
1149: }
1151: /* Show solution */
1152: GVecViewFromOptions(ctx->uExact, "Exact Solution");
1153: GVecViewFromOptions(ctx->u, "Solution");
1154: return(0);
1155: }
1157: #undef __FUNCT__
1159: int StokesComputeFlowField(StokesContext ctx) {
1160: int its; /* The iteration count for the linear solver */
1161: int loop;
1164: /* Get command-line options */
1165: StokesContextSetup(ctx);
1167: for(loop = 0; loop < ctx->numLoops; loop++) {
1168: if (loop >= ctx->refineStep) {
1169: StokesRefineGrid(ctx);
1170: }
1172: /* Setup problem */
1173: StokesSetupGrid(ctx);
1174: StokesCreateStructures(ctx);
1175: StokesSetupStructures(ctx);
1177: /* Check the exact solution */
1178: StokesCheckSolution(ctx, ctx->uExact, ctx->pExact, "exact");
1180: /* Solve system */
1181: StokesSolve(ctx, ctx->f, ctx->u, &its);
1183: /* Check the computed solution */
1184: StokesCheckSolution(ctx, ctx->u, ctx->p, "computed");
1186: /* Cleanup */
1187: StokesDestroyStructures(ctx);
1188: }
1189: return(0);
1190: }