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