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