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