Actual source code: schur.c

  1: #ifdef PETSC_RCS_HEADER
  2: static char vcid[] = "$Id: schur.c,v 1.15 2000/10/04 02:16:21 knepley Exp $";
  3: #endif
  4: /*
  5:   Defines the Schur preconditioner due to Ahmed Sameh for AIJ matrices
  6: */
 7:  #include src/sles/pc/pcimpl.h
 8:  #include schur.h

 10: extern int GridResetConstrainedMultiply_Private(Grid, GMat);

 12: /*
 13:   Math:

 15:   There are several possible options for this type of preconditioner. I call it
 16:   a Schur complement preconditioner because we split the problem into blocks

 18:         /  A   B 
 19:          B^T  0 /

 21:   and generally must provide a mechanism for solving (preconditioning) the
 22:   Schur complement B^T A^{-1} B at each itertion of the overall preconditioner.
 23:   We could include Uzawa type preconditioners, but since they are so simple, I
 24:   have already abstracted them using GMatCreateUzawa() and related functions.
 25:     Here I am concerned with preconditioners for the entire saddle point problem
 26:   above, the simplest being a block diagonal (Jacobi)

 28:         / hat A      0   
 29:             0    -hat S /

 31:   or the slightly more sophisticated block triangular

 33:         / hat A      B   
 34:             0    -hat S /

 36:   culminating in the full block LU factorization

 38:         / hat A  0  / I   hat A^{-1} B  
 39:            B^T   I /  0      -hat S    /

 41:   These all depend on solving the Schur complement S, so we need an effective
 42:   preconditioner for it. Elman, Wathen and Silvester propose the "BFBT"
 43:   preconditioner which essentially projects the problem onto the pressure space
 44:   and then applies the convection diffusion operator. So we have

 46:         (B^T A^{-1} B)^{-1} = (B^T B)^{-1} B^T A B (B^T B)^{-1}
 47:                             = A_p (B^T B)^{-1}

 49:   where in the last step they assume a commutation relation of the form

 51:         A B  = B A_p

 53:   so that the momentum equations on the pressure mesh are equivalent in some
 54:   sense to those on the velocity mesh. In addition, they recommend several cycles
 55:   of multigrid for hat A^{-1}.
 56:     Sameh instead proposes to precondition S with the pressure lapalcian L_p and
 57:   perform fast solves with hat A using the balanced method. I do not understand
 58:   why Wathen, et.al. do not do the full block LU instead of the triangular thing,
 59:   but we can test this. We can afford to do an iterative method for S since
 60:   applications of hat A^{-1} are cheap, but this may be compared with the sparse
 61:   factorizations of Saad.
 62:     Lastly, it must be noted that you can solve a preconditioned system for hat A
 63:   or hat S, or merely apply the preconditioners. The cost of solving must be
 64:   weighed against the extra outer iterations incured by the weaker approximation.

 66:   Programming:

 68:     We allow several modes for this preconditioner since parts of the matrix may
 69:   not be available, e.g. if the matrix is given as implicitly P^T A P. The flag
 70:   explicitConstraints signals that a constrained matrix has been explicitly formed.
 71:   Otherwise, we use GridGetConstraintMatrix() to retrieve the projector P onto the
 72:   constraint space.
 73:     Also, we may have new variables introduced by constraints. These will be present
 74:   already in explicit matrices, but are not yet included in the implicit formulation.
 75:   Thus, we have functions which add this contribution to a given vector, with functions
 76:   like GVecEvaluateJacobianConstrained(). This function operates in the projected space
 77:   so that we get an action like

 79:         P^T A P + G

 81:   where G incorporates the contribution of the new fields.
 82: */

 84: /*--------------------------------------------- Public Functions ----------------------------------------------------*/
 85: #undef  __FUNCT__
 87: /*@C PCSchurSetGradientOperator
 88:         This function sets the operator and its transpose, which define the
 89:   constraint manifold in the saddle-point problem. For instance, in
 90:   the Stokes problem, this operator is the gradient, and its transpose
 91:   is the divergence. It must be called prior to PCSetUp().

 93:   Collective on PC

 95:   Input Parameters:
 96: + pc     - The preconditioner context
 97: . gradOp - The gradient operator
 98: - divOp  - The divergence operator

100:   Level: intermediate

102: .keywords schur, set, gradient
103: .seealso PCSchurGetIterationNumber()
104: @*/
105: int PCSchurSetGradientOperator(PC pc, int gradOp, int divOp)
106: {
107:   PC_Schur *s;

111:   s = (PC_Schur *) pc->data;
112:   s->gradOp = gradOp;
113:   s->divOp  = divOp;
114:   return(0);
115: }
116: #undef  __FUNCT__
118: /*@C
119:   PCSchurGetIterationNumber - Get the iteration numbers since the solve was started.

121:   Not collective

123:   Input Parameter:
124: . pc       - The preconditioner context

126:   Output
127: + its      - The total number of iterations used to solve A
128: - schurIts - The number of iterations used to solve the Schur complement

130:   Level: intermediate

132: .keywords schur, iteration
133: .seealso PCSchurSetGradientOperator()
134: @*/
135: int PCSchurGetIterationNumber(PC pc, int *its, int *schurIts)
136: {
137:   PC_Schur *s;

143:   s = (PC_Schur *) pc->data;
144:   *its      = s->iter;
145:   *schurIts = s->schurIter;
146:   return(0);
147: }

149: #undef __FUNCT__  
151: /*@C
152:   PCSchurInnerMonitor - Print the residual norm of the solver for the inner system A.

154:   Collective on KSP

156:   Input Parameters:
157: . ksp   - The KSP context
158: . n     - The iteration number
159: . rnorm - The 2-norm of residual
160: . ctx   - The Schur context

162:   Level: intermediate

164: .keywords: KSP, schur, monitor, norm
165: .seealso: KSPSetMonitor(), PCSchurMonitor()
166: @*/
167: int PCSchurInnerMonitor(KSP ksp, int n, PetscReal rnorm, void *ctx)
168: {
169:   MPI_Comm comm;
170:   int      ierr;

174:   PetscObjectGetComm((PetscObject) ksp, &comm);
175:   PetscPrintf(comm, "  inner iter = %d, residual norm %g n", n, rnorm);
176:   return(0);
177: }

179: #undef __FUNCT__  
181: /*@C
182:   PCSchurMonitor - Print the residual norm of the solver for the Schur complement system S.

184:   Collective on KSP

186:   Input Parameters:
187: . ksp   - The KSP context
188: . n     - The iteration number
189: . rnorm - The 2-norm of residual
190: . ctx   - The Schur context

192:   Level: intermediate

194: .keywords: KSP, schur, monitor, norm
195: .seealso: KSPSetMonitor(), PCSchurInnerMonitor()
196: @*/
197: int PCSchurMonitor(KSP ksp, int n, PetscReal rnorm, void *ctx)
198: {
199:   MPI_Comm comm;
200:   int      ierr;

204:   PetscObjectGetComm((PetscObject) ksp, &comm);
205:   PetscPrintf(comm, "  schur iter = %d, residual norm %g n", n, rnorm);
206:   return(0);
207: }


210: #undef __FUNCT__  
212: /*@C
213:   PCSchurSolveMonitor - Print the number of iterates for each inner solve in the Schur complement system.

215:   Collective on KSP

217:   Input Parameters:
218: . ksp   - The KSP context
219: . n     - The iteration number
220: . rnorm - The 2-norm of residual
221: . ctx   - The Schur context

223:   Level: intermediate

225: .keywords: KSP, schur, monitor, norm
226: .seealso: KSPSetMonitor(), PCSchurMonitor(), PCSchurInnerMonitor()
227: @*/
228: int PCSchurSolveMonitor(KSP ksp, int n, PetscReal rnorm, void *ctx)
229: {
230:   PC_Schur *s = (PC_Schur *) ctx;
231:   KSP       innerKSP;
232:   int       its;
233:   int       ierr;

237:   SLESGetKSP(s->sles, &innerKSP);
238:   KSPGetIterationNumber(innerKSP, &its);
239:   if (n > 0)
240:     PetscLogInfo(ksp, "PCSchur: %d iterations in A solve %dn", its, n);
241:   return(0);
242: }

244: /*--------------------------------------------- Private Functions ---------------------------------------------------*/
245: #undef __FUNCT__  
247: static int PCDestroyStructures_Schur(PC pc)
248: {
249:   PC_Schur *s = (PC_Schur *) pc->data;
250:   int       ierr;

253:   /* Destroy operator structures */
254:   PetscFree(s->momOps);
255:   PetscFree(s->momOpIsALE);
256:   PetscFree(s->momOpAlphas);

258:   /* Destroy variable orderings */
259:   VarOrderingDestroy(s->sOrder);
260:   VarOrderingDestroy(s->tOrder);
261:   LocalVarOrderingDestroy(s->sLocOrder);
262:   LocalVarOrderingDestroy(s->tLocOrder);

264:   /* Destroy matrices and reorderings */
265:   GMatDestroyUzawa(s->S);
266:   GMatDestroy(s->A);
267:   if (s->sparseA != PETSC_NULL) {
268:     GMatDestroy(s->sparseA);
269:   }
270:   GMatDestroy(s->B);
271:   if (s->lap != PETSC_NULL) {
272:     GMatDestroy(s->lap);
273:   }
274:   if (s->rowPerm != PETSC_NULL) {
275:     ISDestroy(s->rowPerm);
276:     ISDestroy(s->colPerm);
277:   }
278:   /* Destroy field structures */
279:   VecDestroy(s->u);
280:   VecDestroy(s->uNew);
281:   VecDestroy(s->uNew2);
282:   VecDestroy(s->p);
283:   VecDestroy(s->pNew);
284:   VecScatterDestroy(s->uScatter);
285:   VecScatterDestroy(s->pScatter);
286:   if ((s->isConstrained == PETSC_TRUE) && (s->explicitConstraints == PETSC_FALSE)) {
287:     VecDestroy(s->projX);
288:     VecDestroy(s->projY);
289:     VecDestroy(s->projZ);
290:   }

292:   return(0);
293: }

295: #undef  __FUNCT__
297: int PCValidQ_Schur(PC pc)
298: {
299:   PC_Schur *s = (PC_Schur *) pc->data;

302:   if (pc->setupcalled < 2)
303:     PetscFunctionReturn(1);

306:   /* Check dimensions */

308:   /* Check fields */
314:   if ((s->isConstrained == PETSC_TRUE) && (s->explicitConstraints == PETSC_FALSE)) {
318:   }

322:   /* Check inner solvers */

326:   /* Check matrices */
330:   if (s->sparseA != PETSC_NULL) {
332:   }

334:   return(0);
335: }

337: #undef __FUNCT__  
339: static int PCDebug_Schur(PC pc)
340: {
341:   PC_Schur    *s = (PC_Schur *) pc->data;
342:   Vec          x, y;
343:   Mat          P;
344:   Vec          xx;
345:   Vec          yy;
346:   PetscScalar *xArray;
347:   PetscScalar  zero = 0.0, one = 1.0, minusOne = -1.0;
348:   PetscReal    uNorm, pNorm;
349:   int          row, size;
350:   int          ierr;

353:   VecDuplicate(pc->vec, &x);
354:   VecDuplicate(pc->vec, &y);
355:   if (s->explicitConstraints == PETSC_FALSE) {
356:     GridGetConstraintMatrix(s->grid, &P);
357:     xx   = s->projX;
358:     yy   = s->projY;
359:   } else {
360:     xx   = x;
361:     yy   = y;
362:   }

364:   /* Check the matrix-vector product */
365:   VecGetSize(x, &size);
366:   VecGetArray(x, &xArray);
367:   for(row = 0; row < size; row++) {
368:     /* Start with e_row */
369:     VecSet(&zero, x);
370:     xArray[row] = 1.0;
371:     /* hat A x */
372:     MatMult(pc->mat, x, y);
373:     if (s->explicitConstraints == PETSC_FALSE) {
374:       /* Project vectors */
375:       MatMult(P, x, xx);
376:       MatMult(P, y, yy);
377:     }
378:     /* New hat A x */
379:     VecScatterBegin(xx, s->u, INSERT_VALUES, SCATTER_FORWARD, s->uScatter);
380:     VecScatterEnd(xx, s->u, INSERT_VALUES, SCATTER_FORWARD, s->uScatter);
381:     VecScatterBegin(xx, s->p, INSERT_VALUES, SCATTER_FORWARD, s->pScatter);
382:     VecScatterEnd(xx, s->p, INSERT_VALUES, SCATTER_FORWARD, s->pScatter);
383:     if (s->colPerm != PETSC_NULL) {
384:       VecPermute(s->u, s->colPerm, PETSC_FALSE);
385:     }
386:     /* uNew2 = hat A u + B p, pNew = B u */
387:     MatMult(s->A, s->u, s->uNew);
388:     MatMult(s->B, s->p, s->uNew2);
389:     VecAXPY(&one, s->uNew2, s->uNew);
390:     MatMultTranspose(s->B, s->u, s->pNew);
391:     if (s->rowPerm != PETSC_NULL) {
392:       VecPermute(s->uNew, s->rowPerm, PETSC_TRUE);
393:     }
394:     if (s->explicitConstraints == PETSC_FALSE) {
395:       /* Add in contribution from new variables */
396:       /* GVecEvaluateJacobianConstrained(xx, xx);                                                     */
397:       VecScatterBegin(xx, s->uNew, ADD_VALUES, SCATTER_FORWARD, s->uScatter);
398:       VecScatterEnd(xx, s->uNew, ADD_VALUES, SCATTER_FORWARD, s->uScatter);
399:       VecScatterBegin(xx, s->pNew, ADD_VALUES, SCATTER_FORWARD, s->pScatter);
400:       VecScatterEnd(xx, s->pNew, ADD_VALUES, SCATTER_FORWARD, s->pScatter);
401:     }
402:     /* Compare field vectors */
403:     VecScatterBegin(yy, s->u, INSERT_VALUES, SCATTER_FORWARD, s->uScatter);
404:     VecScatterEnd(yy, s->u, INSERT_VALUES, SCATTER_FORWARD, s->uScatter);
405:     VecScatterBegin(yy, s->p, INSERT_VALUES, SCATTER_FORWARD, s->pScatter);
406:     VecScatterEnd(yy, s->p, INSERT_VALUES, SCATTER_FORWARD, s->pScatter);
407:     VecAXPY(&minusOne, s->u, s->uNew);
408:     VecAXPY(&minusOne, s->p, s->pNew);
409:     VecNorm(s->uNew, NORM_2, &uNorm);
410:     VecNorm(s->pNew, NORM_2, &pNorm);
411:     if ((uNorm > 1.0e-8) || (pNorm > 1.0e-8)) {
412:       SETERRQ(PETSC_ERR_PLIB, "Invalid field matrix");
413:     }
414:   }

416:   VecRestoreArray(x, &xArray);
417:   VecDestroy(x);
418:   VecDestroy(y);
419:   return(0);
420: }

422: #undef __FUNCT__  
424: int PCSetFromOptions_Inner(PC pc, SLES sles, SLES schurSles)
425: {
426:   PC_Schur  *s = (PC_Schur *) pc->data;
427:   KSP        ksp;
428:   char      *prefix;
429:   PetscTruth opt;
430:   int        ierr;

433:   SLESGetOptionsPrefix(sles, &prefix);
434:   SLESGetKSP(sles, &ksp);
435:   PetscOptionsHasName(prefix, "-ksp_monitor", &opt);
436:   if (opt == PETSC_TRUE) {
437:     KSPClearMonitor(ksp);
438:     KSPSetMonitor(ksp, PCSchurInnerMonitor, s, PETSC_NULL);
439:   }
440:   SLESGetOptionsPrefix(schurSles, &prefix);
441:   SLESGetKSP(schurSles, &ksp);
442:   PetscOptionsHasName(prefix, "-ksp_monitor", &opt);
443:   if (opt == PETSC_TRUE) {
444:     KSPClearMonitor(ksp);
445:     KSPSetMonitor(ksp, PCSchurMonitor, s, PETSC_NULL);
446:   }
447:   PetscOptionsHasName(prefix, "-ksp_solve_monitor", &opt);
448:   if (opt == PETSC_TRUE) {
449:     KSPClearMonitor(ksp);
450:     KSPSetMonitor(ksp, PCSchurSolveMonitor, s, PETSC_NULL);
451:   }
452:   return(0);
453: }

455: #undef __FUNCT__  
457: static int PCSetFromOptions_Schur(PC pc)
458: {
459:   PC_Schur  *s = (PC_Schur *) pc->data;
460:   PetscTruth opt;
461:   int        ierr;

464:   /* Enable explicit constraints */
465:   PetscOptionsHasName(pc->prefix, "-pc_schur_explicit", &opt);
466:   if (opt == PETSC_TRUE)
467:     s->explicitConstraints = PETSC_TRUE;

469:   /* Enable laplacian preconditioning of the Schur complement */
470:   PetscOptionsHasName(pc->prefix, "-pc_schur_lap", &opt);
471:   if (opt == PETSC_TRUE)
472:     s->useLaplacian = PETSC_TRUE;

474:   /* Enable Mathematica support */
475:   PetscOptionsHasName(pc->prefix, "-pc_schur_math", &opt);
476:   if (opt == PETSC_TRUE) {
477:     s->useMath = PETSC_TRUE;
478:     PetscViewerMathematicaOpen(pc->comm, PETSC_DECIDE, PETSC_NULL, PETSC_NULL, &s->mathViewer);
479:   }

481:   /* Setup inner solvers */
482:   PCSetFromOptions_Inner(pc, s->sles, s->schurSles);
483:   return(0);
484: }

486: #undef __FUNCT__  
488: static int PCSchurGetMomentumOperators(PC pc)
489: {
490:   PC_Schur   *s = (PC_Schur *) pc->data;
491:   int         numOps;         /* Grid op query: The number of grid matrix operators used for A */
492:   int         sField, tField; /* Grid op query: The shape and test function fields */
493:   PetscScalar alpha;          /* Grid op query: The scalar multiplier */
494:   PetscTruth  isALE;          /* Grid op query: The flag for ALE operators */
495:   int         op;
496:   int         ierr;

499:   GridGetNumMatOperators(s->grid, &s->numMomOps);
500:   s->numMomOps  -= 2;
501:   PetscMalloc(s->numMomOps * sizeof(int),        &s->momOps);
502:   PetscMalloc(s->numMomOps * sizeof(PetscTruth), &s->momOpIsALE);
503:   PetscMalloc(s->numMomOps * sizeof(double),     &s->momOpAlphas);
504:   /* Find gradient fields */
505:   ierr   = GridGetMatOperatorStart(s->grid, &op, &sField, &tField, &alpha, &isALE);
506:   numOps = 0;
507:   while(op >= 0) {
508:     if (op == s->gradOp) {
509:       s->sField      = sField;
510:       s->tField      = tField;
511:       s->gradOpAlpha = alpha;
512:     } else if (op != s->divOp) {
513:       s->momOps[numOps]      = op;
514:       s->momOpIsALE[numOps]  = isALE;
515:       s->momOpAlphas[numOps] = alpha;
516:       numOps++;
517:     }
518:     GridGetMatOperatorNext(s->grid, &op, &sField, &tField, &alpha, &isALE);
519:   }
520:   if ((s->sField < 0) || (s->tField < 0)) {
521:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Gradient operator not found in grid.");
522:   }
523:   if (numOps != s->numMomOps) {
524:     SETERRQ(PETSC_ERR_PLIB, "Inconsistent operator setup in grid");
525:   }

527:   return(0);
528: }

530: #undef __FUNCT__  
532: static int PCSchurReorderMatrices(PC pc)
533: {
534:   PC_Schur  *s = (PC_Schur *) pc->data;
535:   Mat        newA, newB;     /* The reordering matrices */
536:   char       orderType[256]; /* The type of reordering */
537:   IS         tempIS;         /* The identity ordering for B */
538:   int        bw;             /* The bandwidth limitation on the reordered matrix */
539:   PetscReal  frac;           /* The fractional bandwidth limitation on the reordered matrix */
540:   double     tol;            /* The drop tolerance */
541:   int        m, n;
542:   PetscTruth opt, issparse;
543:   int        ierr;

546:   /* Reduce the bandwidth of A */
547:   PetscOptionsHasName(pc->prefix, "-pc_schur_reorder", &opt);
548:   if (opt == PETSC_TRUE) {
549:     MatGetOrdering(s->A, MATORDERING_RCM, &s->rowPerm, &s->colPerm);
550:     PetscOptionsGetString(pc->prefix, "-pc_schur_reorder", orderType, 255, &opt);
551:     if (opt == PETSC_TRUE) {
552:       PetscStrcasecmp(orderType, "sparse", &issparse);
553:       if (issparse) {
554:         bw   = PETSC_DECIDE;
555:         PetscOptionsGetInt(pc->prefix, "-pc_schur_sparsify_bw", &bw, &opt);
556:         frac = 0.0;
557:         PetscOptionsGetReal(pc->prefix, "-pc_schur_sparsify_frac", &frac, &opt);
558:         tol  = 0.0;
559:         PetscOptionsGetReal(pc->prefix, "-pc_schur_sparsify_tol", &frac, &opt);
560:         MatPermuteSparsify(s->A, bw, frac, tol, s->rowPerm, s->colPerm, &s->sparseA);
561:       }
562:       MatPermute(s->A, s->rowPerm, s->colPerm, &newA);
563:     }
564:     MatDestroy(s->A);
565:     MatGetSize(s->B, &m, &n);
566:     ISCreateStride(pc->comm, n, 0, 1, &tempIS);
567:     ISSetPermutation(tempIS);
568:     MatPermute(s->B, s->rowPerm, tempIS, &newB);
569:     ISDestroy(tempIS);
570:     MatDestroy(s->B);
571:     s->A = newA;
572:     s->B = newB;
573:     PetscObjectCompose((PetscObject) s->A, "Grid", (PetscObject) s->grid);
574:     PetscObjectCompose((PetscObject) s->B, "Grid", (PetscObject) s->grid);
575:   }

577:   /* View the new matrices */
578:   PetscOptionsHasName(PETSC_NULL, "-pc_schur_view", &opt);
579:   if (opt == PETSC_TRUE) {
580:     PetscViewer viewer = PETSC_VIEWER_DRAW_(pc->comm);
581:     PetscDraw   draw;

583:     PetscViewerDrawGetDraw(viewer, 0, &draw);
584:     PetscDrawSetPause(draw, -1);
585:     MatView(s->A, viewer);
586:     if (s->sparseA != PETSC_NULL) {
587:       MatView(s->sparseA, viewer);
588:     }
589:     MatView(s->B, viewer);
590:   }

592:   return(0);
593: }

595: #undef __FUNCT__  
597: static int PCSchurCreateMatrices(PC pc)
598: {
599:   PC_Schur     *s    = (PC_Schur *) pc->data;
600:   Grid          grid = s->grid;
601:   VarOrdering   tempOrder;
602:   int           op;
603:   int           ierr;

606:   GridIsConstrained(grid, &s->isConstrained);
607:   /* Create test function variable orderings */
608:   VarOrderingCreateGeneral(grid, 1, &s->tField, &s->tOrder);
609:   LocalVarOrderingCreate(grid, 1, &s->tField, &s->tLocOrder);

611:   /* Create matrices */
612:   if ((s->isConstrained == PETSC_FALSE) || (s->explicitConstraints == PETSC_FALSE)) {
613:     /* Create shape function variable orderings */
614:     VarOrderingCreateGeneral(grid, 1, &s->sField, &s->sOrder);
615:     LocalVarOrderingCreate(grid, 1, &s->sField, &s->sLocOrder);
616:     GMatCreateRectangular(grid, s->tOrder, s->tOrder, &s->A);
617:     GMatCreateRectangular(grid, s->sOrder, s->tOrder, &s->B);
618:     if (s->useLaplacian == PETSC_TRUE) {
619:       GMatCreateRectangular(grid, s->sOrder, s->sOrder, &s->lap);
620:     }
621:   } else {
622:     /* Constrain test function variable ordering */
623:     VarOrderingConstrain(grid, s->tOrder, &tempOrder);
624:     VarOrderingDestroy(s->tOrder);
625:     s->tOrder = tempOrder;
626:     /* Create shape function variable orderings */
627:     VarOrderingCreateGeneral(grid, 1, &s->sField, &s->sOrder);
628:     LocalVarOrderingCreate(grid, 1, &s->sField, &s->sLocOrder);
629:     GMatCreateRectangular(grid, s->tOrder, s->tOrder, &s->A);
630:     GMatCreateRectangular(grid, s->sOrder, s->tOrder, &s->B);
631:     if (s->useLaplacian == PETSC_TRUE) {
632:       GMatCreateRectangular(grid, s->sOrder, s->sOrder, &s->lap);
633:     }
634:   }
635:   MatSetOption(s->A, MAT_NEW_NONZERO_ALLOCATION_ERR);
636:   MatSetOption(s->B, MAT_NEW_NONZERO_ALLOCATION_ERR);
637:   if (s->useLaplacian == PETSC_TRUE) {
638:     MatSetOption(s->lap, MAT_NEW_NONZERO_ALLOCATION_ERR);
639:   }

641:   /* Construct the operators */
642:   if ((s->isConstrained == PETSC_FALSE) || (s->explicitConstraints == PETSC_FALSE)) {
643:     for(op = 0; op < s->numMomOps; op++) {
644:       if (s->momOpIsALE[op] == PETSC_FALSE) {
645:         GMatEvaluateOperatorGalerkin(s->A, PETSC_NULL, 1, &s->tField, &s->tField, s->momOps[op], s->momOpAlphas[op],
646:                                             MAT_FINAL_ASSEMBLY, s);
647: 
648:       } else {
649:         GMatEvaluateALEOperatorGalerkin(s->A, 1, &s->tField, s->tOrder, s->tLocOrder, &s->tField, s->tOrder,
650:                                                s->tLocOrder, s->momOps[op], s->momOpAlphas[op],
651:                                                MAT_FINAL_ASSEMBLY, s);
652: 
653:       }
654:     }
655:     GMatEvaluateOperatorGalerkin(s->B, PETSC_NULL, 1, &s->sField, &s->tField, s->gradOp, s->gradOpAlpha, MAT_FINAL_ASSEMBLY, s);
656: 
657:     if (s->useLaplacian == PETSC_TRUE) {
658:       GMatEvaluateOperatorGalerkin(s->lap, PETSC_NULL, 1, &s->sField, &s->sField, LAPLACIAN, 1.0, MAT_FINAL_ASSEMBLY, s);
659: 
660:     }
661:   } else {
662:     for(op = 0; op < s->numMomOps; op++) {
663:       if (s->momOpIsALE[op] == PETSC_FALSE) {
664:         GMatEvaluateOperatorGalerkin(s->A, PETSC_NULL, 1, &s->tField, &s->tField, s->momOps[op], s->momOpAlphas[op],
665:                                             MAT_FINAL_ASSEMBLY, s);
666: 
667:       } else {
668:         GMatEvaluateALEConstrainedOperatorGalerkin(s->A, 1, &s->tField, s->tOrder, s->tLocOrder, &s->tField, s->tOrder,
669:                                                           s->tLocOrder, s->momOps[op], s->momOpAlphas[op], MAT_FINAL_ASSEMBLY, s);
670: 
671:       }
672:     }
673:     GMatEvaluateNewFields(s->A, 1, &s->tField, s->tOrder, s->tLocOrder, &s->tField, s->tOrder, s->tLocOrder, 1.0,
674:                                  MAT_FINAL_ASSEMBLY);
675: 
676:     GMatEvaluateOperatorGalerkin(s->B, PETSC_NULL, 1, &s->sField, &s->tField, s->gradOp, s->gradOpAlpha, MAT_FINAL_ASSEMBLY, s);
677: 
678:     if (s->useLaplacian == PETSC_TRUE) {
679:       GMatEvaluateOperatorGalerkin(s->lap, PETSC_NULL, 1, &s->sField, &s->sField, LAPLACIAN, 1.0, MAT_FINAL_ASSEMBLY, s);
680: 
681:     }
682:     /* Reset matrix multiplication */
683:     GridResetConstrainedMultiply_Private(grid, s->A);
684:     GridResetConstrainedMultiply_Private(grid, s->B);
685:     if (s->useLaplacian == PETSC_TRUE) {
686:       GridResetConstrainedMultiply_Private(grid, s->lap);
687:     }
688:   }

690:   /* Reduce the bandwidth of A */
691:   PCSchurReorderMatrices(pc);

693:   return(0);
694: }

696: #undef __FUNCT__  
698: static int PCSchurCreateInnerSLES(PC pc)
699: {
700:   PC_Schur *s = (PC_Schur *) pc->data;
701:   int       ierr;

704:   /* Create the momentum solver */
705:   SLESCreate(pc->comm, &s->sles);
706:   SLESAppendOptionsPrefix(s->sles, "pc_schur_inner_");

708:   /* Create the Schur complement solver */
709:   SLESCreate(pc->comm, &s->schurSles);
710:   SLESAppendOptionsPrefix(s->schurSles, "pc_schur_");
711:   return(0);
712: }

714: #undef __FUNCT__  
716: static int PCSchurSetupInnerSLES(PC pc)
717: {
718:   PC_Schur *s = (PC_Schur *) pc->data;
719:   KSP       innerKsp;
720:   PC        innerPc;
721:   int       ierr;

724:   /* Create the momentum solver */
725:   SLESCreate(pc->comm, &s->sles);
726:   if (s->sparseA != PETSC_NULL) {
727:     SLESSetOperators(s->sles, s->A, s->sparseA, DIFFERENT_NONZERO_PATTERN);
728:   } else {
729:     SLESSetOperators(s->sles, s->A, s->A,       SAME_NONZERO_PATTERN);
730:   }
731:   SLESAppendOptionsPrefix(s->sles, "pc_schur_inner_");
732:   SLESGetKSP(s->sles, &innerKsp);
733:   SLESGetPC(s->sles, &innerPc);
734:   KSPSetType(innerKsp, KSPPREONLY);
735:   PCSetType(innerPc, PCLU);
736:   SLESSetFromOptions(s->sles);

738:   /* Create the Schur complement solver */
739:   SLESCreate(pc->comm, &s->schurSles);
740:   GMatCreateUzawa(s->A, s->B, s->sles, &s->S);
741:   SLESSetOperators(s->schurSles, s->S, s->S, SAME_NONZERO_PATTERN);
742:   SLESAppendOptionsPrefix(s->schurSles, "pc_schur_");
743:   SLESGetKSP(s->schurSles, &innerKsp);
744:   SLESGetPC(s->schurSles, &innerPc);
745:   KSPSetType(innerKsp, KSPGMRES);
746:   if (s->useLaplacian == PETSC_TRUE) {
747:     PCSetOperators(innerPc, s->S, s->lap, DIFFERENT_NONZERO_PATTERN);
748:     PCSetType(innerPc, PCLU);
749:   } else {
750:     PCSetType(innerPc, PCNONE);
751:   }
752:   SLESSetFromOptions(s->schurSles);

754:   /* Setup monitors */
755:   PCSetFromOptions_Inner(pc, s->sles, s->schurSles);
756:   return(0);
757: }

759: #undef __FUNCT__  
761: static int PCCreateStructures_Schur(PC pc)
762: {
763:   PC_Schur *s    = (PC_Schur *) pc->data;
764:   Grid      grid = s->grid;
765:   int       ierr;

768:   /* Create matrices */
769:   PCSchurCreateMatrices(pc);

771:   /* Create the work vectors and field scatters */
772:   GVecCreateRectangular(grid, s->tOrder, &s->u);
773:   GVecDuplicate(s->u, &s->uNew);
774:   GVecDuplicate(s->u, &s->uNew2);
775:   GVecCreateRectangular(grid, s->sOrder, &s->p);
776:   GVecDuplicate(s->p, &s->pNew);
777:   if ((s->isConstrained == PETSC_TRUE) && (s->explicitConstraints == PETSC_FALSE)) {
778:     GVecCreate(grid, &s->projX);
779:     GVecDuplicate(s->projX, &s->projY);
780:     GVecCreateConstrained(grid, &s->projZ);
781:     GVecScatterCreate(s->projX, s->u, &s->uScatter);
782:     GVecScatterCreate(s->projX, s->p, &s->pScatter);
783:   } else {
784:     GVecScatterCreate(pc->vec, s->u, &s->uScatter);
785:     GVecScatterCreate(pc->vec, s->p, &s->pScatter);
786:   }

788:   /* Setup Inner Solvers */
789:   PCSchurSetupInnerSLES(pc);

791:   return(0);
792: }

794: #undef __FUNCT__  
796: static int PCSetUp_Schur(PC pc)
797: {
798:   PC_Schur  *s = (PC_Schur *) pc->data;
799:   PetscTruth opt;
800:   int        ierr;

803:   /* This lets the factorization become stale, that is uses the same preconditioner for several Newton steps */
804:   if (pc->setupcalled > 0) {
805:     return(0);
806:   }
807:   pc->setupcalled = 2;

809:   /* Get the grid */
810:   GVecGetGrid(pc->vec, &s->grid);

812:   /* Divide up the problem */
813:   PCSchurGetMomentumOperators(pc);

815:   /* Create matrices, vector storage, scatters, and inner solvers */
816:   PCCreateStructures_Schur(pc);

818: #ifdef PETSC_USE_BOPT_g
819:   PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);
820: #endif
821:   PetscOptionsHasName(PETSC_NULL, "-pc_schur_debug", &opt);
822:   if (opt == PETSC_TRUE) {
823:     PCDebug_Schur(pc);
824:   }

826:   return(0);
827: }

829: #undef __FUNCT__  
831: static int PCPreSolve_Schur(PC pc, KSP ksp, Vec x,  Vec rhs)
832: {
833:   PC_Schur *s = (PC_Schur *) pc->data;

836:   s->iter      = 0;
837:   s->schurIter = 0;
838:   return(0);
839: }

841: #undef __FUNCT__  
843: static int PCPostSolve_Schur(PC pc, KSP ksp, Vec x, Vec rhs)
844: {
846:   return(0);
847: }

849: #undef __FUNCT__  
851: static int PCApply_Schur(PC pc, Vec x, Vec y)
852: {
853:   PC_Schur   *s  = (PC_Schur *) pc->data;
854:   Vec         xx = x;
855:   Vec         yy = y;
856:   Vec         z  = s->projZ;
857:   Vec         tempSol;
858:   Mat         P, invP;
859:   PetscScalar minusOne = -1.0;
860:   int         its, schurIts;
861:   PetscTruth  opt;
862:   int         ierr;

865:   PetscOptionsHasName(PETSC_NULL, "-pc_schur_test", &opt);
866:   if (opt == PETSC_TRUE) {
867:     /* Set the input vector to Ax */
868:     VecDuplicate(x, &tempSol);
869:     VecCopy(x, tempSol);
870:     MatMult(pc->mat, x, y);
871:     VecCopy(y, x);
872:   }

874:   /* Project: Apply (P^+)^T = P (P^T P)^{-1} */
875:   if ((s->isConstrained == PETSC_TRUE) && (s->explicitConstraints == PETSC_FALSE)) {
876:     xx   = s->projX;
877:     yy   = s->projY;
878:     GridGetConstraintMatrix(s->grid, &P);
879:     GridGetConstraintInverse(s->grid, &invP);
880:     MatMult(invP, x, z);
881:     MatMult(P, z, xx);
882:   }

884:   /*            /    A^{-1}    0      0    / u    /     A^{-1} u        / w_u 
885:      L^{-1} x = | -B^T A^{-1}  I      0   | | p | = | p - B^T A^{-1} u | = | w_p |
886:                       0       0   G^{-1} /  U /        G^{-1} U     /    w_U /
887:   */
888:   VecScatterBegin(xx, s->u, INSERT_VALUES, SCATTER_FORWARD, s->uScatter);
889:   VecScatterEnd(xx, s->u, INSERT_VALUES, SCATTER_FORWARD, s->uScatter);
890:   /* Permute u */
891:   if (s->colPerm != PETSC_NULL) {
892:     VecPermute(s->u, s->colPerm, PETSC_FALSE);
893:   }
894:   /* Velocity solve: Now w_u = s->uNew */
895:   SLESSolve(s->sles, s->u, s->uNew, &its);
896:   s->iter += its;
897:   PetscLogInfo(pc, "PCSchur: %d iterations in A solve %dn", its, 0);
898:   /* Pressure solve: Now w_p = s->p */
899:   MatMultTranspose(s->B, s->uNew, s->p);
900:   VecScale(&minusOne, s->p);
901:   VecScatterBegin(xx, s->p, ADD_VALUES, SCATTER_FORWARD, s->pScatter);
902:   VecScatterEnd(xx, s->p, ADD_VALUES, SCATTER_FORWARD, s->pScatter);

904:   /*            /   I    A^{-1} B (B^T A^{-1} B)^{-1}   0    / w_u    / w_u + A^{-1} B (B^T A^{-1} B)^{-1} w_p 
905:      U^{-1} w = |   0       -(B^T A^{-1} B)^{-1}        0   | | w_p | = |         -(B^T A^{-1} B)^{-1} w_p       |
906:                    0                 0                 I   /  w_U /                         w_U               /
907:   */
908:   /* Schur solve: Now s->pNew = -S^{-1} w_p */
909:   SLESSolve(s->schurSles, s->p, s->pNew, &schurIts);
910:   VecScale(&minusOne, s->pNew);
911:   s->schurIter += schurIts;
912:   PetscLogInfo(pc, "PCSchur: %d iterations in B^T A^{-1} B solven", schurIts);
913:   /* Velocity solve: Now s->uNew2 = -A^{-1} B S^{-1} w_p */
914:   MatMult(s->B, s->pNew, s->u);
915:   SLESSolve(s->sles, s->u, s->uNew2, &its);
916:   s->iter += its;
917:   PetscLogInfo(pc, "PCSchur: %d iterations in A solve %dn", its, schurIts+1);
918:   /* Now s->uNew = w_u - A^{-1} B S^{-1} w_p */
919:   VecAXPY(&minusOne, s->uNew2, s->uNew);
920:   /* Permute uNew */
921:   if (s->rowPerm != PETSC_NULL) {
922:     VecPermute(s->uNew, s->rowPerm, PETSC_TRUE);
923:   }
924:   /* Put u and p in y */
925:   VecScatterBegin(s->uNew, yy, INSERT_VALUES, SCATTER_REVERSE, s->uScatter);
926:   VecScatterEnd(s->uNew, yy, INSERT_VALUES, SCATTER_REVERSE, s->uScatter);
927:   VecScatterBegin(s->pNew, yy, INSERT_VALUES, SCATTER_REVERSE, s->pScatter);
928:   VecScatterEnd(s->pNew, yy, INSERT_VALUES, SCATTER_REVERSE, s->pScatter);

930:   /* Project back: Apply P^+ = (P^T P)^{-1} P^T */
931:   if ((s->isConstrained == PETSC_TRUE) && (s->explicitConstraints == PETSC_FALSE)) {
932:     /* Project solution back onto constrained unknowns */
933:     MatMultTranspose(P, yy, z);
934:     /* Particle solve: Add G^{-1} U = w_U to z */
935:     GVecSolveJacobianConstrained(z, x);
936:     /* Projection Solve: Now y = (P^T P)^{-1} z */
937:     MatMult(invP, z, y);
938:   }

940:   PetscOptionsHasName(PETSC_NULL, "-pc_schur_test", &opt);
941:   if (opt == PETSC_TRUE) {
942:     PetscReal solNorm;

944:     /* Check that we get A^{-1} A x = x */
945:     VecAXPY(&minusOne, y, tempSol);
946:     VecNorm(tempSol, NORM_2, &solNorm);
947:     VecDestroy(tempSol);
948:     PetscPrintf(pc->comm, "Error in test solution: %gn", solNorm);
949:   }
950:   return(0);
951: }

953: #undef __FUNCT__  
955: static int PCApplyTrans_Schur(PC pc, Vec x, Vec y)
956: {
958:   return(0);
959: }

961: #undef __FUNCT__  
963: static int PCApplySymmetricLeft_Schur(PC pc, Vec x, Vec y)
964: {
966:   return(0);
967: }

969: #undef __FUNCT__  
971: static int PCApplySymmetricRight_Schur(PC pc, Vec x, Vec y)
972: {
974:   return(0);
975: }

977: #undef __FUNCT__  
979: static int PCDestroy_Schur(PC pc)
980: {
981:   PC_Schur *s = (PC_Schur *) pc->data;
982:   int       ierr;

985:   /* Destroy inner solvers */
986:   SLESDestroy(s->sles);
987:   SLESDestroy(s->schurSles);
988:   if (pc->setupcalled) {
989:     PCDestroyStructures_Schur(pc);
990:   }
991:   if (s->mathViewer != PETSC_NULL) {
992:     PetscViewerDestroy(s->mathViewer);
993:   }
994:   PetscFree(s);
995:   return(0);
996: }

998: EXTERN_C_BEGIN
999: #undef __FUNCT__  
1001: int PCCreate_Schur(PC pc)
1002: {
1003:   PC_Schur *s;
1004:   int       ierr;

1007:   PetscNew(PC_Schur, &s);
1008:   PetscLogObjectMemory(pc, sizeof(PC_Schur));
1009:   pc->data                = (void *) s;

1011:   pc->ops->setup               = PCSetUp_Schur;
1012:   pc->ops->apply               = PCApply_Schur;
1013:   pc->ops->applyrichardson     = PETSC_NULL;
1014:   pc->ops->applyBA             = PETSC_NULL;
1015:   pc->ops->applytranspose      = PCApplyTrans_Schur;
1016:   pc->ops->applyBAtranspose    = PETSC_NULL;
1017:   pc->ops->setfromoptions      = PCSetFromOptions_Schur;
1018:   pc->ops->presolve            = PCPreSolve_Schur;
1019:   pc->ops->postsolve           = PCPostSolve_Schur;
1020:   pc->ops->getfactoredmatrix   = PETSC_NULL;
1021:   pc->ops->applysymmetricleft  = PCApplySymmetricLeft_Schur;
1022:   pc->ops->applysymmetricright = PCApplySymmetricRight_Schur;
1023:   pc->ops->setuponblocks       = PETSC_NULL;
1024:   pc->ops->destroy             = PCDestroy_Schur;
1025:   pc->ops->view                = PETSC_NULL;
1026:   pc->modifysubmatrices        = PETSC_NULL;

1028:   s->grid                 = PETSC_NULL;
1029:   s->isConstrained        = PETSC_FALSE;
1030:   s->explicitConstraints  = PETSC_FALSE;
1031:   s->useMath              = PETSC_FALSE;
1032:   s->mathViewer           = PETSC_NULL;

1034:   s->u                    = PETSC_NULL;
1035:   s->uNew                 = PETSC_NULL;
1036:   s->uNew2                = PETSC_NULL;
1037:   s->p                    = PETSC_NULL;
1038:   s->pNew                 = PETSC_NULL;
1039:   s->uScatter             = PETSC_NULL;
1040:   s->pScatter             = PETSC_NULL;
1041:   s->projX                = PETSC_NULL;
1042:   s->projY                = PETSC_NULL;
1043:   s->projZ                = PETSC_NULL;

1045:   s->sles                 = PETSC_NULL;
1046:   s->schurSles            = PETSC_NULL;
1047:   s->iter                 = 0;
1048:   s->schurIter            = 0;

1050:   s->numMomOps            = -1;
1051:   s->gradOp               = GRADIENT;
1052:   s->divOp                = DIVERGENCE;
1053:   s->gradOpAlpha          = 0.0;
1054:   s->momOps               = PETSC_NULL;
1055:   s->momOpIsALE           = PETSC_NULL;
1056:   s->momOpAlphas          = PETSC_NULL;
1057:   s->sField               = -1;
1058:   s->tField               = -1;
1059:   s->sOrder               = PETSC_NULL;
1060:   s->sLocOrder            = PETSC_NULL;
1061:   s->tOrder               = PETSC_NULL;
1062:   s->tLocOrder            = PETSC_NULL;
1063:   s->A                    = PETSC_NULL;
1064:   s->sparseA              = PETSC_NULL;
1065:   s->B                    = PETSC_NULL;
1066:   s->S                    = PETSC_NULL;
1067:   s->lap                  = PETSC_NULL;
1068:   s->rowPerm              = PETSC_NULL;
1069:   s->colPerm              = PETSC_NULL;

1071:   /* Create the inner solvers */
1072:   PCSchurCreateInnerSLES(pc);

1074:   return(0);
1075: }
1076: EXTERN_C_END