Actual source code: control.c
1: #ifdef PETSC_RCS_HEADER
2: static char vcid[] = "$Id: control.c,v 1.5 2000/02/01 17:02:51 knepley Exp $";
3: #endif
5: static char help[] = "This is the Control problem with Dirchlet boundary conditions.\n\n";
7: int CONTROL_COOKIE;
8: int CONTROL_ComputeField;
10: #include control.h
12: /*
13: Here we are solving the Control equation,
15: L = min 1/2 \int^2_0 u(x)^2 dx subject to \nabla u = c(x), u(0) = 1, and \abs{c} \le 1
16: u
18: In order to enforce the equality constraints, we introduce additional terms in the Lagrangian,
20: L = min 1/2 \int^2_0 dx u^2 + \int^2_0 dx \lambda(x) (\nabla u - c)
21: u
23: and we require that the first variation vanish.
25: / \int^2_0 dx (u - \nabla\lambda) \
26: \delta L = | \int^2_0 dx (\nabla u - c) |
27: \ \int^2_0 dx -\lambda /
29: In order to drive the nonlinear solver, we also need the second variation
31: / I -\nabla 0 \
32: | |
33: | \nabla\cdot 0 -I |
34: | |
35: \ 0 -I 0 /
37: Thus we now let
39: CostFunctional() --> L
40: GradientFunction() --> \delta L
41: */
42: #undef __FUNCT__
44: int main(int argc, char **args) {
45: ControlContext ctx; /* Holds problem specific information */
46: int ierr;
49: PetscInitialize(&argc, &args, 0, help); CHKERRABORT(PETSC_COMM_WORLD, ierr);
50: TaoInitialize(&argc, &args, 0, help); CHKERRABORT(PETSC_COMM_WORLD, ierr);
52: ControlContextCreate(PETSC_COMM_WORLD, &ctx); CHKERRABORT(PETSC_COMM_WORLD, ierr);
53: ControlComputeField(ctx); CHKERRABORT(PETSC_COMM_WORLD, ierr);
54: ControlContextDestroy(ctx); CHKERRABORT(PETSC_COMM_WORLD, ierr);
56: CHKMEMQ;
57: PetscFinalize();
58: TaoFinalize();
59: return(0);
60: }
64: /*@C
65: ControlInitializePackage - This function initializes everything in the Control package.
67: Input Parameter:
68: path - The dynamic library path, or PETSC_NULL
70: Level: developer
72: .keywords: Control, initialize, package
73: .seealso: PetscInitialize()
74: @*/
75: int ControlInitializePackage(char *path) {
76: static PetscTruth initialized = PETSC_FALSE;
77: char logList[256];
78: char *className;
79: PetscTruth opt;
80: int ierr;
83: if (initialized == PETSC_TRUE) return(0);
84: initialized = PETSC_TRUE;
85: /* Register Classes */
86: PetscLogClassRegister(&CONTROL_COOKIE, "Control");
87: /* Register Constructors and Serializers */
88: /* Register Events */
89: PetscLogEventRegister(&CONTROL_ComputeField, "ControlComputeField", CONTROL_COOKIE);
90: /* Process info exclusions */
91: PetscOptionsGetString(PETSC_NULL, "-log_info_exclude", logList, 256, &opt);
92: if (opt == PETSC_TRUE) {
93: PetscStrstr(logList, "ControlContext", &className);
94: if (className) {
95: PetscLogInfoDeactivateClass(CONTROL_COOKIE);
96: }
97: }
98: /* Process summary exclusions */
99: PetscOptionsGetString(PETSC_NULL, "-log_summary_exclude", logList, 256, &opt);
100: if (opt == PETSC_TRUE) {
101: PetscStrstr(logList, "ControlContext", &className);
102: if (className) {
103: PetscLogEventDeactivateClass(CONTROL_COOKIE);
104: }
105: }
106: return(0);
107: }
109: /*-------------------------------------------- ControlContext Creation ------------------------------------------------*/
110: #undef __FUNCT__
112: /*@
113: ControlContextCreate - This function initializes the Control context.
115: Collective on MPI_Comm
117: Input Parameter:
118: . comm - The communicator
120: Output Parameter:
121: . sCtx - The ControlContext
123: Level: beginner
125: .keywords: Control, context, create
126: .seealso: ControlContextDestroy(), ControlContextPrint(), ControlContextSetup()
127: @*/
128: int ControlContextCreate(MPI_Comm comm, ControlContext *sCtx) {
129: ControlContext ctx;
132: /* Setup context */
133: PetscHeaderCreate(ctx, _ControlContext, int, PETSC_VIEWER_COOKIE, -1, "ControlContext", comm, 0, 0);
134: PetscLogObjectCreate(ctx);
135: PetscLogObjectMemory(ctx, sizeof(struct _ControlContext));
137: /* Initialize subobjects */
138: ctx->grid = PETSC_NULL;
139: ctx->sles = PETSC_NULL;
140: ctx->A = PETSC_NULL;
141: ctx->u = PETSC_NULL;
142: ctx->f = PETSC_NULL;
143: ctx->uExact = PETSC_NULL;
144: /* Setup domain */
145: ctx->geometryCtx.size[0] = 2.0;
146: ctx->geometryCtx.size[1] = 2.0;
147: ctx->geometryCtx.start[0] = 0.0;
148: ctx->geometryCtx.start[1] = 0.0;
149: /* Setup refinement */
150: ctx->geometryCtx.maxArea = 0.5;
151: ctx->geometryCtx.areaFunc = PointFunctionConstant;
152: ctx->geometryCtx.areaCtx = PETSC_NULL;
153: /* Initialize problem loop */
154: ctx->dim = 2;
155: ctx->linear = PETSC_FALSE;
156: ctx->refineStep = 0;
157: ctx->numLoops = 0;
159: *sCtx = ctx;
160: return(0);
161: }
163: #undef __FUNCT__
165: /*@
166: ControlContextDestroy - This function destroys the Control context.
168: Collective on ControlContext
170: Input Parameter:
171: . ctx - The ControlContext
173: Level: beginner
175: .keywords: Control, context, destroy
176: .seealso: ControlContextCreate(), ControlContextSetup()
177: @*/
178: int ControlContextDestroy(ControlContext ctx) {
179: Grid grid = ctx->grid;
180: int ierr;
183: if (--ctx->refct > 0) SETERRQ(PETSC_ERR_PLIB, "Control context should not be referenced more than once");
184: TaoDestroy(ctx->tao);
185: TaoApplicationDestroy(ctx->taoApp);
186: GridFinalizeBC(grid);
187: GridDestroy(grid);
188: PetscLogObjectDestroy(ctx);
189: PetscHeaderDestroy(ctx);
190: return(0);
191: }
193: /*--------------------------------------------- ControlContext Setup -------------------------------------------------*/
194: #undef __FUNCT__
196: int ControlContextSetup(ControlContext ctx) {
197: PetscTruth opt;
198: int ierr;
201: /* Determine the problem dimension */
202: PetscOptionsGetInt(PETSC_NULL, "-dim", &ctx->dim, &opt);
203: /* Determine the element type */
204: PetscOptionsHasName(PETSC_NULL, "-linear", &ctx->linear);
205: /* Determine how many systems to solve */
206: PetscOptionsGetInt(PETSC_NULL, "-num_systems", &ctx->numLoops, &opt);
207: /* The first iteration at which to refine the mesh */
208: PetscOptionsGetInt(PETSC_NULL, "-refine_step", &ctx->refineStep, &opt);
209: /* The maximum area of any triangle in the refined mesh */
210: PetscOptionsGetReal("mesh", "-max_area", &ctx->geometryCtx.maxArea, &opt);
212: /* Create main problem */
213: ControlContextCreateGrid(ctx);
214: ControlContextCreateTao(ctx);
215: return(0);
216: }
218: /*------------------------------------------------ Grid Creation -----------------------------------------------------*/
219: #undef __FUNCT__
221: /*@
222: ControlContextCreateMeshBoundary - This function creates a mesh boundary for the main problem.
224: Collective on ControlContext
226: Input Parameter:
227: . ctx - A ControlContext with problem specific information
229: Level: beginner
231: .seealso ControlContextDestroyMeshBoundary(), ControlContextCreateMesh()
232: @*/
233: int ControlContextCreateMeshBoundary(ControlContext ctx) {
234: MPI_Comm comm;
235: MeshGeometryContext *geomCtx = &ctx->geometryCtx;
236: int ierr;
239: PetscObjectGetComm((PetscObject) ctx, &comm);
240: switch(ctx->dim) {
241: case 1:
242: MeshBoundary1DCreateSimple(comm, geomCtx, &ctx->boundaryCtx);
243: break;
244: case 2:
245: MeshBoundary2DCreateSimple(comm, geomCtx, &ctx->boundaryCtx);
246: break;
247: default:
248: SETERRQ1(PETSC_ERR_SUP, "Do not support meshes of dimension %d", ctx->dim);
249: }
250: return(0);
251: }
253: #undef __FUNCT__
255: /*@
256: ControlContextDestroyMeshBoundary - This function destroys the mesh boundary for the main problem.
258: Collective on ControlContext
260: Input Parameter:
261: . ctx - A ControlContext with problem specific information
263: Level: beginner
265: .seealso ControlContextCreateMeshBoundary(), ControlContextCreateMesh()
266: @*/
267: int ControlContextDestroyMeshBoundary(ControlContext ctx) {
268: MPI_Comm comm;
269: int ierr;
272: PetscObjectGetComm((PetscObject) ctx, &comm);
273: switch(ctx->dim) {
274: case 1:
275: MeshBoundary1DDestroy(comm, &ctx->boundaryCtx);
276: break;
277: case 2:
278: MeshBoundary2DDestroy(comm, &ctx->boundaryCtx);
279: break;
280: default:
281: SETERRQ1(PETSC_ERR_SUP, "Do not support meshes of dimension %d", ctx->dim);
282: }
283: return(0);
284: }
286: #undef __FUNCT__
288: /*@
289: ControlContextCreateMesh - This function creates a mesh for the main problem.
291: Collective on ControlContext
293: Input Parameter:
294: . ctx - A ControlContext with problem specific information
296: Output Parameter:
297: . m - The Mesh
299: Options Database Keys:
300: . mesh_refine - Refines the mesh based on area criteria
301: . mesh_max_area - The maximum area of an element
303: Level: beginner
305: .seealso ControlRefineGrid(), ControlDestroyGrid()
306: @*/
307: int ControlContextCreateMesh(ControlContext ctx, Mesh *m) {
308: MeshGeometryContext *geomCtx = &ctx->geometryCtx;
309: MPI_Comm comm;
310: Mesh mesh;
311: Partition part;
312: int totalElements, totalNodes, totalEdges;
313: char name[1024];
314: int d;
315: int ierr;
318: PetscObjectGetComm((PetscObject) ctx, &comm);
319: ControlContextCreateMeshBoundary(ctx);
320: MeshCreate(comm, &mesh);
321: MeshSetDimension(mesh, ctx->dim);
322: for(d = 0; d < ctx->dim; d++) {
323: MeshSetPeriodicDimension(mesh, d, geomCtx->isPeriodic[d]);
324: }
325: switch(ctx->dim) {
326: case 1:
327: MeshSetNumCorners(mesh, 3);
328: break;
329: case 2:
330: if (ctx->linear == PETSC_TRUE) {
331: MeshSetNumCorners(mesh, 4);
332: } else {
333: MeshSetNumCorners(mesh, 6);
334: }
335: break;
336: default:
337: SETERRQ1(PETSC_ERR_SUP, "Do not support meshes of dimension %d", ctx->dim);
338: }
339: MeshSetBoundary(mesh, &ctx->boundaryCtx);
340: sprintf(name, "mesh.r%.4g", geomCtx->maxArea);
341: PetscObjectSetName((PetscObject) mesh, name);
342: MeshSetFromOptions(mesh);
343: ControlContextDestroyMeshBoundary(ctx);
344: /* Setup refinement */
345: if (ctx->geometryCtx.areaCtx != PETSC_NULL) {
346: MeshSetUserContext(mesh, geomCtx->areaCtx);
347: } else {
348: MeshSetUserContext(mesh, &geomCtx->maxArea);
349: }
350: /* Report on mesh */
351: MeshGetPartition(mesh, &part);
352: switch(ctx->dim) {
353: case 1:
354: PartitionGetTotalElements(part, &totalElements);
355: PartitionGetTotalNodes(part, &totalNodes);
356: PetscPrintf(comm, "Elements: %d Nodes: %d\n", totalElements, totalNodes);
357: break;
358: case 2:
359: PartitionGetTotalElements(part, &totalElements);
360: PartitionGetTotalNodes(part, &totalNodes);
361: PartitionGetTotalEdges(part, &totalEdges);
362: PetscPrintf(comm, "Elements: %d Nodes: %d Edges: %d\n", totalElements, totalNodes, totalEdges);
363: break;
364: default:
365: SETERRQ1(PETSC_ERR_SUP, "Do not support meshes of dimension %d", ctx->dim);
366: }
368: *m = mesh;
369: return(0);
370: }
372: #undef __FUNCT__
374: /*@
375: ControlContextCreateGrid - This function creates a grid for the main problem.
377: Collective on ControlContext
379: Input Parameter:
380: . ctx - A ControlContext with problem specific information
382: Level: beginner
384: .seealso ControlRefineGrid(), ControlDestroyGrid()
385: @*/
386: int ControlContextCreateGrid(ControlContext ctx) {
387: Mesh mesh;
388: Grid grid;
389: int ierr;
392: /* Construct the mesh */
393: ControlContextCreateMesh(ctx, &mesh);
394: /* Construct the grid */
395: GridCreate(mesh, &grid);
396: MeshDestroy(mesh);
397: switch(ctx->dim) {
398: case 1:
399: if (ctx->linear == PETSC_TRUE) {
400: GridAddField(grid, "u", DISCRETIZATION_TRIANGULAR_1D_LINEAR, 1, PETSC_NULL);
401: GridAddField(grid, "lambda", DISCRETIZATION_TRIANGULAR_1D_CONSTANT, 1, PETSC_NULL);
402: GridAddField(grid, "c", DISCRETIZATION_TRIANGULAR_1D_CONSTANT, 1, PETSC_NULL);
403: } else {
404: GridAddField(grid, "u", DISCRETIZATION_TRIANGULAR_1D_QUADRATIC, 1, PETSC_NULL);
405: GridAddField(grid, "lambda", DISCRETIZATION_TRIANGULAR_1D_LINEAR, 1, PETSC_NULL);
406: GridAddField(grid, "c", DISCRETIZATION_TRIANGULAR_1D_LINEAR, 1, PETSC_NULL);
407: }
408: break;
409: case 2:
410: if (ctx->linear == PETSC_TRUE) {
411: GridAddField(grid, "u", DISCRETIZATION_TRIANGULAR_2D_LINEAR, 2, PETSC_NULL);
412: GridAddField(grid, "c", DISCRETIZATION_TRIANGULAR_2D_CONSTANT, 2, PETSC_NULL);
413: } else {
414: GridAddField(grid, "u", DISCRETIZATION_TRIANGULAR_2D_QUADRATIC, 2, PETSC_NULL);
415: GridAddField(grid, "c", DISCRETIZATION_TRIANGULAR_2D_LINEAR, 2, PETSC_NULL);
416: }
417: break;
418: default:
419: SETERRQ1(PETSC_ERR_SUP, "Do not support meshes of dimension %d", ctx->dim);
420: }
421: GridSetFromOptions(grid);
423: ctx->grid = grid;
424: return(0);
425: }
427: #undef __FUNCT__
429: /*@
430: ControlContextCreateTao - This function creates a Tao Application for the main problem.
432: Collective on ControlContext
434: Input Parameter:
435: . ctx - A ControlContext with problem specific information
437: Level: beginner
439: .seealso ControlCreateGrid(), ControlDestroyTao()
440: @*/
441: int ControlContextCreateTao(ControlContext ctx) {
442: TaoMethod method = "tao_blmvm"; /* minimization method */
443: TAO_SOLVER tao; /* TAO_SOLVER solver context */
444: TAO_APPLICATION taoApp;
445: int ierr;
448: TaoCreate(ctx->comm, method, &tao);
449: TaoPetscApplicationCreate(ctx->comm, &taoApp);
451: ctx->tao = tao;
452: ctx->taoApp = taoApp;
453: return(0);
454: }
456: #undef __FUNCT__
458: /*@
459: ControlRefineGrid - This function refines the mesh for the main grid.
461: Collective on ControlContext
463: Input Parameters:
464: + ctx - A ControlContext
466: Options Database Keys:
467: . mesh_max_area - The maximum area na element may have
469: Level: beginner
471: .seealso ControlCreateGRid(), ControlDestroyGrid()
472: @*/
473: int ControlRefineGrid(ControlContext ctx) {
474: Grid oldGrid = ctx->grid;
475: Grid grid;
476: Mesh mesh;
477: Partition part;
478: MeshGeometryContext *geomCtx = &ctx->geometryCtx;
479: int totalElements, totalNodes, totalEdges;
480: int ierr;
483: /* Construct a refined mesh */
484: if (geomCtx->areaCtx != PETSC_NULL) {
485: GridRefineMesh(oldGrid, geomCtx->areaFunc, &grid);
486: GridGetMesh(grid, &mesh);
487: MeshSetUserContext(mesh, geomCtx->areaCtx);
488: } else {
489: geomCtx->maxArea *= 0.5;
490: GridRefineMesh(oldGrid, geomCtx->areaFunc, &grid);
491: GridGetMesh(grid, &mesh);
492: MeshSetUserContext(mesh, &geomCtx->maxArea);
493: }
494: GridSetOptionsPrefix(grid, "ref_");
495: GridSetFromOptions(grid);
497: MeshGetPartition(mesh, &part);
498: PartitionGetTotalElements(part, &totalElements);
499: PartitionGetTotalNodes(part, &totalNodes);
500: PartitionGetTotalEdges(part, &totalEdges);
501: PetscPrintf(ctx->comm, "Elements: %d Nodes: %d Edges: %d\n", totalElements, totalNodes, totalEdges);
502: CHKMEMQ;
504: /* Replace old grid with refined grid */
505: GridFinalizeBC(oldGrid);
506: GridDestroy(oldGrid);
507: ctx->grid = grid;
508: return(0);
509: }
511: /*------------------------------------------------- Grid Setup -------------------------------------------------------*/
512: #undef __FUNCT__
514: /*@
515: ControlSetupGrid - This function sets all the functions,
516: operators , and boundary conditions for the problem.
517: It also sets the parameters associated with the fields.
519: Collective on Grid
521: Input Parameters:
522: . grid - The grid
523: . ctx - A ControlContext
525: Options Database Keys:
526: . use_laplacian - Use the Laplacian in stead of the Rate-of-Strain tensor
528: Level: intermediate
530: .seealso ControlCreateGrid()
531: @*/
532: int ControlSetupGrid(ControlContext ctx) {
533: Grid grid = ctx->grid;
534: int uField = 0;
535: int lField = 1;
536: int cField = 2;
537: int ierr;
540: /* Setup Problem */
541: GridSetFieldName(grid, uField, "State");
542: GridSetFieldName(grid, lField, "Multiplier");
543: GridSetFieldName(grid, cField, "Control");
544: GridSetActiveField(grid, uField);
545: GridAddActiveField(grid, lField);
546: GridAddActiveField(grid, cField);
548: /* Setup Gradient */
549: GridSetRhsOperator(grid, uField, uField, IDENTITY, 1.0, PETSC_FALSE, PETSC_NULL);
550: GridAddRhsOperator(grid, lField, uField, GRADIENT, -1.0, PETSC_FALSE, PETSC_NULL);
551: GridAddRhsOperator(grid, uField, lField, DIVERGENCE, 1.0, PETSC_FALSE, PETSC_NULL);
552: GridAddRhsOperator(grid, cField, lField, IDENTITY, -1.0, PETSC_FALSE, PETSC_NULL);
553: GridAddRhsOperator(grid, lField, cField, IDENTITY, -1.0, PETSC_FALSE, PETSC_NULL);
554: ControlSetupRhsFunction(grid, ctx);
556: /* Setup Hessian */
557: GridSetMatOperator(grid, uField, uField, IDENTITY, 1.0, PETSC_FALSE, PETSC_NULL);
558: GridAddMatOperator(grid, lField, uField, GRADIENT, -1.0, PETSC_FALSE, PETSC_NULL);
559: GridAddMatOperator(grid, uField, lField, DIVERGENCE, 1.0, PETSC_FALSE, PETSC_NULL);
560: GridAddMatOperator(grid, cField, lField, IDENTITY, -1.0, PETSC_FALSE, PETSC_NULL);
561: GridAddMatOperator(grid, lField, cField, IDENTITY, -1.0, PETSC_FALSE, PETSC_NULL);
563: /* Setup Dirchlet boundary conditions */
564: ControlSetupBC(grid, ctx);
565: return(0);
566: }
568: #undef __FUNCT__
570: /*@ ControlSetupRhsFunction
571: ControlSetupRhsFunction - This function chooses a forcing function for the problem.
573: Collective on Grid
575: Input Parameters:
576: . grid - The grid
577: . ctx - A ControlContext
579: Level: intermediate
581: Options Database Keys:
583: .seealso ControlSetupGrid
584: @*/
585: int ControlSetupRhsFunction(Grid grid, ControlContext ctx) {
587: #if 0
588: GridAddRhsFunction(grid, uField, RhsFunction, 1.0);
589: #endif
590: return(0);
591: }
593: #undef __FUNCT__
595: /*@ ControlSetupBC
596: ControlSetupBC - This function chooses boundary conditions for the problem.
598: Collective on Grid
600: Input Parameters:
601: . grid - The grid
602: . ctx - A ControlContext
604: Level: intermediate
606: Options Database Keys:
607: . bc_reduce - Explicitly reduce the system using boundary conditions
609: .seealso ControlSetupGrid()
610: @*/
611: int ControlSetupBC(Grid grid, ControlContext ctx) {
612: int uField = 0;
613: PetscTruth reduceSystem, reduceElement;
614: int ierr;
617: PetscOptionsHasName(PETSC_NULL, "-bc_reduce", &reduceSystem);
618: PetscOptionsHasName(PETSC_NULL, "-bc_reduce_elem", &reduceElement);
619: switch(ctx->dim) {
620: case 1:
621: GridAddBC(grid, BOTTOM_BD, uField, PointFunctionOne, reduceSystem);
622: break;
623: case 2:
624: break;
625: default:
626: SETERRQ1(PETSC_ERR_SUP, "Do not support meshes of dimension %d", ctx->dim);
627: }
628: GridSetReduceElement(grid, reduceElement);
629: GridSetBCContext(grid, ctx);
630: return(0);
631: }
633: #undef __FUNCT__
635: /*@
636: ControlSetupTao - This function sets the function, gradient, and hessian evaluators, as well as bounds.
638: Collective on ControlContext
640: Input Parameter:
641: . ctx - A ControlContext
643: Level: intermediate
645: .seealso ControlCreateTao()
646: @*/
647: int ControlSetupTao(ControlContext ctx) {
651: ControlSetupBounds(ctx);
652: GVecDuplicate(ctx->u, &ctx->lagrangianDensity);
654: TaoSetPetscFunctionGradient(ctx->taoApp, ctx->u, ctx->f, ControlFormFunctionGradient, ctx);
655: TaoSetPetscHessian(ctx->taoApp, ctx->A, ctx->A, ControlFormHessian, ctx);
656: TaoSetPetscVariableBounds(ctx->taoApp, ctx->lowerBound, ctx->upperBound);
657: TaoSetApplication(ctx->tao, ctx->taoApp);
659: TaoSetFromOptions(ctx->tao);
660: return(0);
661: }
663: #undef __FUNCT__
665: /*@
666: ControlSetupBounds - This function sets the variable bounds.
668: Collective on ControlContext
670: Input Parameter:
671: . ctx - A ControlContext
673: Level: intermediate
675: .seealso ControlSetupTao()
676: @*/
677: int ControlSetupBounds(ControlContext ctx) {
678: int uField = 0;
679: int lField = 1;
680: int cField = 2;
684: VecDuplicate(ctx->u, &ctx->lowerBound);
685: VecDuplicate(ctx->u, &ctx->upperBound);
686: #if 0
687: GVecEvaluateFunction(ctx->lowerBound, 1, &uField, PointFunctionOne, PETSC_MIN, PETSC_NULL);
688: GVecEvaluateFunction(ctx->lowerBound, 1, &lField, PointFunctionOne, PETSC_MIN, PETSC_NULL);
689: GVecEvaluateFunction(ctx->lowerBound, 1, &cField, PointFunctionZero, 1.0, PETSC_NULL);
690: GVecEvaluateFunction(ctx->upperBound, 1, &uField, PointFunctionOne, PETSC_MAX, PETSC_NULL);
691: GVecEvaluateFunction(ctx->upperBound, 1, &lField, PointFunctionOne, PETSC_MAX, PETSC_NULL);
692: GVecEvaluateFunction(ctx->upperBound, 1, &cField, PointFunctionOne, 1.0, PETSC_NULL);
693: #else
694: GVecEvaluateFunction(ctx->lowerBound, 1, &uField, PointFunctionOne, -1.0e10, PETSC_NULL);
695: GVecEvaluateFunction(ctx->lowerBound, 1, &lField, PointFunctionOne, -1.0e10, PETSC_NULL);
696: GVecEvaluateFunction(ctx->lowerBound, 1, &cField, PointFunctionZero, 1.0, PETSC_NULL);
697: GVecEvaluateFunction(ctx->upperBound, 1, &uField, PointFunctionOne, 1.0e10, PETSC_NULL);
698: GVecEvaluateFunction(ctx->upperBound, 1, &lField, PointFunctionOne, 1.0e10, PETSC_NULL);
699: GVecEvaluateFunction(ctx->upperBound, 1, &cField, PointFunctionOne, 1.0, PETSC_NULL);
700: #endif
701: return(0);
702: }
704: /*------------------------------------------------- SLES Setup -------------------------------------------------------*/
705: #undef __FUNCT__
707: int ControlCreateStructures(ControlContext ctx) {
711: /* Create the linear solver */
712: SLESCreate(ctx->comm, &ctx->sles);
713: SLESSetFromOptions(ctx->sles);
714: /* Create solution, rhs, and exact solution vectors */
715: GVecCreate(ctx->grid, &ctx->u);
716: PetscObjectSetName((PetscObject) ctx->u, "Solution");
717: GVecDuplicate(ctx->u, &ctx->f);
718: PetscObjectSetName((PetscObject) ctx->f, "Rhs");
719: GVecDuplicate(ctx->u, &ctx->uExact);
720: PetscObjectSetName((PetscObject) ctx->uExact, "ExactSolution");
721: /* Create the system matrix */
722: GMatCreate(ctx->grid, &ctx->A);
723: return(0);
724: }
726: #undef __FUNCT__
728: int ControlSetupKSP(KSP ksp, ControlContext ctx) {
729: GVecErrorKSPMonitorCtx *monCtx = &ctx->monitorCtx;
730: PetscViewer v;
731: PetscDraw draw;
732: PetscTruth opt;
733: int ierr;
736: /* Setup convergence monitors */
737: PetscOptionsHasName(PETSC_NULL, "-error_viewer", &opt);
738: if (opt == PETSC_TRUE) {
739: v = PETSC_VIEWER_DRAW_(ctx->comm);
740: PetscViewerSetFormat(v, PETSC_VIEWER_DRAW_LG);
741: PetscViewerDrawGetDraw(v, 0, &draw);
742: PetscDrawSetTitle(draw, "Error");
743: } else {
744: v = PETSC_VIEWER_STDOUT_(ctx->comm);
745: }
746: monCtx->error_viewer = v;
747: monCtx->solution = ctx->uExact;
748: monCtx->norm_error_viewer = PETSC_VIEWER_STDOUT_(ctx->comm);
749: KSPSetMonitor(ksp, GVecErrorKSPMonitor, monCtx, PETSC_NULL);
750: return(0);
751: }
754: #undef __FUNCT__
756: int ControlSetupPC(PC pc, ControlContext ctx) {
757: MPI_Comm comm;
758: PetscScalar alpha;
759: int lField = 1;
760: int ierr;
763: PetscObjectGetComm((PetscObject) pc, &comm);
764: GVecCreate(ctx->grid, &ctx->constantL);
765: GVecEvaluateFunction(ctx->constantL, 1, &lField, PointFunctionOne, 1.0, PETSC_NULL);
766: VecDot(ctx->constantL, ctx->constantL, &alpha);
767: alpha = 1.0/PetscSqrtScalar(alpha);
768: VecScale(&alpha, ctx->constantL);
769: #if 0
770: MatNullSpaceCreate(comm, PETSC_FALSE, 1, &ctx->constantL, &ctx->nullSpace);
771: MatNullSpaceTest(ctx->nullSpace, ctx->A);
772: PCNullSpaceAttach(pc, ctx->nullSpace);
773: MatNullSpaceRemove(ctx->nullSpace, ctx->uExact, PETSC_NULL);
774: MatNullSpaceDestroy(ctx->nullSpace);
775: #endif
776: return(0);
777: }
779: #undef __FUNCT__
781: int ControlSetupStructures(ControlContext ctx) {
782: KSP ksp;
783: PC pc;
784: int uField = 0;
785: int lField = 1;
786: int cField = 2;
787: MatStructure flag;
788: PetscTruth opt;
789: int ierr;
792: /* Evaluate the rhs */
793: GridEvaluateRhs(ctx->grid, PETSC_NULL, ctx->f, (PetscObject) ctx);
794: /* Evaluate the exact solution */
795: GVecEvaluateFunction(ctx->uExact, 1, &uField, SolutionFunction, 1.0, ctx);
796: GVecEvaluateFunction(ctx->uExact, 1, &lField, MultiplierSolutionFunction, 1.0, ctx);
797: GVecEvaluateFunction(ctx->uExact, 1, &cField, ControlSolutionFunction, 1.0, ctx);
798: /* Evaluate the system matrix */
799: flag = DIFFERENT_NONZERO_PATTERN;
800: GridEvaluateSystemMatrix(ctx->grid, PETSC_NULL, &ctx->A, &ctx->A, &flag, (PetscObject) ctx);
801: MatCheckSymmetry(ctx->A);
802: /* Apply Dirchlet boundary conditions */
803: GMatSetBoundary(ctx->A, 1.0, ctx);
804: GVecSetBoundary(ctx->f, ctx);
805: /* Setup the linear solver */
806: SLESSetOperators(ctx->sles, ctx->A, ctx->A, SAME_NONZERO_PATTERN);
807: SLESGetKSP(ctx->sles, &ksp);
808: ControlSetupKSP(ksp, ctx);
809: SLESGetPC(ctx->sles, &pc);
810: ControlSetupPC(pc, ctx);
811: /* View structures */
812: PetscOptionsHasName(PETSC_NULL, "-mat_view", &opt);
813: if (opt == PETSC_TRUE) {
814: MatView(ctx->A, PETSC_VIEWER_STDOUT_(ctx->comm));
815: MatView(ctx->A, PETSC_VIEWER_DRAW_(ctx->comm));
816: }
817: PetscOptionsHasName(PETSC_NULL, "-vec_view", &opt);
818: if (opt == PETSC_TRUE) {
819: VecView(ctx->f, PETSC_VIEWER_STDOUT_(ctx->comm));
820: }
821: return(0);
822: }
824: #undef __FUNCT__
826: int ControlDestroyStructures(ControlContext ctx) {
830: SLESDestroy(ctx->sles);
831: GVecDestroy(ctx->f);
832: GVecDestroy(ctx->u);
833: GVecDestroy(ctx->uExact);
834: GMatDestroy(ctx->A);
835: GVecDestroy(ctx->constantL);
836: GVecDestroy(ctx->lagrangianDensity);
837: GVecDestroy(ctx->lowerBound);
838: GVecDestroy(ctx->upperBound);
839: return(0);
840: }
842: /*----------------------------------------------- Sanity Checks ------------------------------------------------------*/
843: #undef __FUNCT__
845: int MatCheckSymmetry(Mat A) {
846: Mat trA;
847: PetscTruth isSym;
848: int ierr;
851: MatTranspose(A, &trA);
852: MatEqual(A, trA, &isSym);
853: MatDestroy(trA);
854: if (isSym == PETSC_FALSE) PetscFunctionReturn(1);
855: return(0);
856: }
858: #undef __FUNCT__
860: int ControlCheckSolution(ControlContext ctx, GVec u, const char type[]) {
861: GVec r;
862: PetscScalar minusOne = -1.0;
863: PetscScalar norm;
864: int ierr;
867: GVecDuplicate(u, &r);
868: /* A u */
869: MatMult(ctx->A, u, r);
870: /* f - A u^* */
871: VecAYPX(&minusOne, ctx->f, r);
872: /* || f - A u || */
873: VecNorm(r, NORM_2, &norm);
874: PetscPrintf(ctx->comm, "Residual of the %s solution: %g\n", type, norm);
876: GVecDestroy(r);
877: return(0);
878: }
880: /*----------------------------------------------- Problem Callbacks --------------------------------------------------*/
881: #undef __FUNCT__
883: /*@
884: SolutionFunction - This function is the velocity solution function for the problem.
886: Not collective
888: Input Parameters:
889: + n - The number of points
890: . comp - The number of components
891: . x,y,z - The points
892: . values - The output
893: - ctx - A ControlContext
895: Level: beginner
897: Note:
898: $ The solution for one dimension is u = 1 - x if x < 1
899: $ 0 if x >= 1
900: $ The solution for two dimensions is u = .
902: .keywords velocity, solution
903: .seealso RhsFunction()
904: @*/
905: int SolutionFunction(int n, int comp, double *x, double *y, double *z, PetscScalar *values, void *ctx) {
906: int i;
909: switch(comp) {
910: case 1:
911: for(i = 0; i < n; i++) {
912: if (x[i] < 1.0) {
913: values[i] = 1.0 - x[i];
914: } else {
915: values[i] = 0.0;
916: }
917: }
918: break;
919: case 2:
920: for(i = 0; i < n; i++) {
921: values[i*2+0] = x[i];
922: values[i*2+1] = y[i];
923: }
924: break;
925: default:
926: SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid number of components %d", comp);
927: }
928: return(0);
929: }
931: #undef __FUNCT__
933: /*@
934: MultiplierSolutionFunction - This function is the solution function for the multiplier \lambda.
936: Not collective
938: Input Parameters:
939: + n - The number of points
940: . comp - The number of components
941: . x,y,z - The points
942: . values - The output
943: - ctx - A ControlContext
945: Level: beginner
947: Note:
948: $ The solution for one dimension is \lambda = 1 if x < 1
949: 0 if x >= 1
950: $ The solution for two dimensions is \lambda =
952: .keywords velocity, solution
953: .seealso RhsFunction()
954: @*/
955: int MultiplierSolutionFunction(int n, int comp, double *x, double *y, double *z, PetscScalar *values, void *ctx) {
956: int i;
959: switch(comp) {
960: case 1:
961: for(i = 0; i < n; i++) {
962: if (x[i] < 1.0) {
963: values[i] = 1.0;
964: } else {
965: values[i] = 0.0;
966: }
967: }
968: break;
969: case 2:
970: for(i = 0; i < n; i++) {
971: values[i*2+0] = x[i];
972: values[i*2+1] = y[i];
973: }
974: break;
975: default:
976: SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid number of components %d", comp);
977: }
978: return(0);
979: }
981: #undef __FUNCT__
983: /*@
984: ControlSolutionFunction - This function is the solution function for the control c.
986: Not collective
988: Input Parameters:
989: + n - The number of points
990: . comp - The number of components
991: . x,y,z - The points
992: . values - The output
993: - ctx - A ControlContext
995: Level: beginner
997: Note:
998: $ The solution for one dimension is c = -1 if x < 1
999: 0 if x >= 1
1000: $ The solution for two dimensions is c =
1002: .keywords velocity, solution
1003: .seealso RhsFunction()
1004: @*/
1005: int ControlSolutionFunction(int n, int comp, double *x, double *y, double *z, PetscScalar *values, void *ctx) {
1006: int i;
1009: switch(comp) {
1010: case 1:
1011: for(i = 0; i < n; i++) {
1012: if (x[i] < 1.0) {
1013: values[i] = -1.0;
1014: } else {
1015: values[i] = 0.0;
1016: }
1017: }
1018: break;
1019: case 2:
1020: for(i = 0; i < n; i++) {
1021: values[i*2+0] = x[i];
1022: values[i*2+1] = y[i];
1023: }
1024: break;
1025: default:
1026: SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid number of components %d", comp);
1027: }
1028: return(0);
1029: }
1031: #undef __FUNCT__
1033: /*@
1034: RhsFunction - This function is the forcing function for the problem.
1036: Not collective
1038: Input Parameters:
1039: + n - The number of points
1040: . comp - The number of components
1041: . x,y,z - The points
1042: . values - The output
1043: - ctx - The ControlContext
1045: Level: beginner
1047: Note:
1048: $ The rhs for one dimension is .
1049: $ The rhs for two dimensions is .
1051: .keywords velocity, rhs
1052: .seealso SolutionFunction()
1053: @*/
1054: int RhsFunction(int n, int comp, double *x, double *y, double *z, PetscScalar *values, void *ctx) {
1055: int i;
1058: switch(comp) {
1059: case 1:
1060: for(i = 0; i < n; i++) {
1061: values[i] = -2.0;
1062: }
1063: break;
1064: case 2:
1065: for(i = 0; i < n; i++) {
1066: values[i*2+0] = -2.0;
1067: values[i*2+1] = 0.0;
1068: }
1069: break;
1070: default:
1071: SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid number of components %d", comp);
1072: }
1073: return(0);
1074: }
1078: /*@
1079: ControlLagrangian - This function returns the Lagrangian L
1081: Not collective
1083: Input Parameters:
1084: + n - The number of points
1085: . comp - The number of components
1086: . x,y,z - The points
1087: . numArgs - The number of inputs
1088: . fieldVals - The field values for each input
1089: - ctx - The ControlContext
1091: Output Parameter:
1092: . values - The output
1094: Level: developer
1096: Note:
1097: $ L for one dimension is 1/2 \int^2_0 dx u^2 - \lambda (\dot x - c).
1098: $ L for two dimensions is .
1100: .keywords Lagrangian
1101: .seealso SolutionFunction()
1102: @*/
1103: int ControlLagrangian(int n, int comp, double *x, double *y, double *z, int numArgs, PetscScalar **fieldVals, PetscScalar *values, void *ctx) {
1104: #if 0
1105: ControlContext cCtx = (ControlContext) ctx;
1106: int dim = cCtx->dim;
1107: int c;
1108: #endif
1109: int i;
1112: if (comp != 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Only one component is supported");
1113: for (i = 0; i < n; i++) {
1114: values[i] = 0.5*fieldVals[0][0]*fieldVals[0][0];
1115: #if 0
1116: for(c = 0; c < comp; c++) {
1117: values[i] = 0.5*fieldVals[0][c*(dim+1)]*fieldVals[0][c*(dim+1)] + \lambda*(fieldVals[0][c*(dim+1)+(c+1)] - control);
1118: }
1119: #endif
1120: }
1121: return(0);
1122: }
1126: /*@
1127: ControlFormFunctionGradient - Evaluates function and corresponding gradient.
1129: Collective on ControlContext
1131: Input Parameters:
1132: + tao - The TAO_SOLVER context
1133: . x - The input vector
1134: - ctx - The ControlContext
1135:
1136: Output Parameters:
1137: + L - The function value
1138: - g - The gradient vector
1140: Level: intermediate
1142: .keywords function, gradient
1143: .seealso ControlFormHessian()
1144: @*/
1145: int ControlFormFunctionGradient(TAO_SOLVER tao, Vec x, double *L, Vec g, void *ctx) {
1146: ControlContext cCtx = (ControlContext) ctx;
1147: GVec l = cCtx->lagrangianDensity;
1148: int fields[3] = {0, 1, 2};
1149: int ierr;
1152: GVecSetBoundary(x, ctx);
1153: GridEvaluateRhs(cCtx->grid, x, g, (PetscObject) ctx);
1154: #if 0
1155: GVecEvaluateNonlinearOperatorGalerkin(l, 1, &x, 3, fields, ControlLagrangian, 1.0, PETSC_FALSE, ctx);
1156: #else
1157: GVecEvaluateNonlinearOperatorGalerkin(l, 1, &x, 1, fields, ControlLagrangian, 1.0, PETSC_FALSE, ctx);
1158: #endif
1159: VecSum(l, L);
1160: return(0);
1161: }
1165: /*@
1166: ControlFormHessian - Evaluates the Hessian matrix.
1168: Collective on ControlContext
1170: Input Parameters:
1171: + tao - The TAO_SOLVER context
1172: . x - The input vector
1173: - ctx - The ControlContext
1175: Output Parameters:
1176: + H - The Hessian matrix
1177: . Hpre - The preconditioning matrix
1178: - flag - The flag indicating matrix structure
1180: Level: intermediate
1182: .keywords hessian
1183: .seealso ControlFormFunctionGradient()
1184: @*/
1185: int ControlFormHessian(TAO_SOLVER tao, Vec x, Mat *H, Mat *Hpre, MatStructure *flag, void *ctx) {
1186: ControlContext cCtx = (ControlContext) ctx;
1187: int ierr;
1190: GVecSetBoundary(x, ctx);
1191: GridEvaluateSystemMatrix(cCtx->grid, x, H, Hpre, flag, (PetscObject) ctx);
1192: MatCheckSymmetry(*H);
1193: /* Indicate that this matrix has the same sparsity pattern during successive iterations;
1194: setting this flag can save significant work in computing the preconditioner for some methods. */
1195: *flag = SAME_NONZERO_PATTERN;
1196: return(0);
1197: }
1199: /*----------------------------------------------- Main Computation ---------------------------------------------------*/
1200: #undef __FUNCT__
1202: int ControlSolve(ControlContext ctx, GVec f, GVec u, int *its) {
1203: PetscReal L; /* The Lagrangian value */
1204: PetscReal gnorm; /* The gradient norm (or duality gap) */
1205: PetscReal cnorm; /* ???Some distance from the constraint manifold */
1206: TaoTerminateReason reason;
1207: int iter;
1208: int ierr;
1211: #if 0
1212: /* Solve A u = f */
1213: SLESSolve(ctx->sles, f, u, its);
1214: #else
1215: /* Solve the bound constrained optimization problem */
1216: TaoSolve(ctx->tao);
1217: TaoGetIterationData(ctx->tao, &iter, &L, &gnorm, &cnorm, 0, &reason);
1218: PetscPrintf(ctx->comm, "Tao Iteration: %d, f: %4.2e, Residual: %4.2e, Infeas: %4.2e\n", iter, L, gnorm, cnorm);
1219: #endif
1221: /* Show solution */
1222: GVecViewFromOptions(ctx->uExact, "Exact Solution");
1223: GVecViewFromOptions(ctx->u, "Solution");
1224: return(0);
1225: }
1227: #undef __FUNCT__
1229: int ControlComputeField(ControlContext ctx) {
1230: int its; /* The iteration count for the linear solver */
1231: int loop;
1234: /* Get command-line options */
1235: ControlContextSetup(ctx);
1237: for(loop = 0; loop < ctx->numLoops; loop++) {
1238: if (loop >= ctx->refineStep) {
1239: ControlRefineGrid(ctx);
1240: }
1242: /* Setup problem */
1243: ControlSetupGrid(ctx);
1244: ControlCreateStructures(ctx);
1245: ControlSetupTao(ctx);
1246: ControlSetupStructures(ctx);
1248: /* Check the exact solution */
1249: ControlCheckSolution(ctx, ctx->uExact, "exact");
1251: /* Solve system */
1252: ControlSolve(ctx, ctx->f, ctx->u, &its);
1254: /* Check the computed solution */
1255: ControlCheckSolution(ctx, ctx->u, "computed");
1257: /* Cleanup */
1258: ControlDestroyStructures(ctx);
1259: }
1260: return(0);
1261: }