Actual source code: lsc.c

  1: #include <petsc/private/pcimpl.h>

  3: typedef struct {
  4:   PetscBool allocated;
  5:   PetscBool scalediag;
  6:   KSP       kspL;
  7:   Vec       scale;
  8:   Vec       x0, y0, x1;
  9:   Mat       L; /* keep a copy to reuse when obtained with L = A10*A01 */
 10: } PC_LSC;

 12: static PetscErrorCode PCLSCAllocate_Private(PC pc)
 13: {
 14:   PC_LSC *lsc = (PC_LSC *)pc->data;
 15:   Mat     A;

 17:   if (lsc->allocated) return 0;
 18:   KSPCreate(PetscObjectComm((PetscObject)pc), &lsc->kspL);
 19:   KSPSetErrorIfNotConverged(lsc->kspL, pc->erroriffailure);
 20:   PetscObjectIncrementTabLevel((PetscObject)lsc->kspL, (PetscObject)pc, 1);
 21:   KSPSetType(lsc->kspL, KSPPREONLY);
 22:   KSPSetOptionsPrefix(lsc->kspL, ((PetscObject)pc)->prefix);
 23:   KSPAppendOptionsPrefix(lsc->kspL, "lsc_");
 24:   MatSchurComplementGetSubMatrices(pc->mat, &A, NULL, NULL, NULL, NULL);
 25:   MatCreateVecs(A, &lsc->x0, &lsc->y0);
 26:   MatCreateVecs(pc->pmat, &lsc->x1, NULL);
 27:   if (lsc->scalediag) VecDuplicate(lsc->x0, &lsc->scale);
 28:   lsc->allocated = PETSC_TRUE;
 29:   return 0;
 30: }

 32: static PetscErrorCode PCSetUp_LSC(PC pc)
 33: {
 34:   PC_LSC *lsc = (PC_LSC *)pc->data;
 35:   Mat     L, Lp, B, C;

 37:   PCLSCAllocate_Private(pc);
 38:   PetscObjectQuery((PetscObject)pc->mat, "LSC_L", (PetscObject *)&L);
 39:   if (!L) PetscObjectQuery((PetscObject)pc->pmat, "LSC_L", (PetscObject *)&L);
 40:   PetscObjectQuery((PetscObject)pc->pmat, "LSC_Lp", (PetscObject *)&Lp);
 41:   if (!Lp) PetscObjectQuery((PetscObject)pc->mat, "LSC_Lp", (PetscObject *)&Lp);
 42:   if (!L) {
 43:     MatSchurComplementGetSubMatrices(pc->mat, NULL, NULL, &B, &C, NULL);
 44:     if (!lsc->L) {
 45:       MatMatMult(C, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &lsc->L);
 46:     } else {
 47:       MatMatMult(C, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &lsc->L);
 48:     }
 49:     Lp = L = lsc->L;
 50:   }
 51:   if (lsc->scale) {
 52:     Mat Ap;
 53:     MatSchurComplementGetSubMatrices(pc->mat, NULL, &Ap, NULL, NULL, NULL);
 54:     MatGetDiagonal(Ap, lsc->scale); /* Should be the mass matrix, but we don't have plumbing for that yet */
 55:     VecReciprocal(lsc->scale);
 56:   }
 57:   KSPSetOperators(lsc->kspL, L, Lp);
 58:   KSPSetFromOptions(lsc->kspL);
 59:   return 0;
 60: }

 62: static PetscErrorCode PCApply_LSC(PC pc, Vec x, Vec y)
 63: {
 64:   PC_LSC *lsc = (PC_LSC *)pc->data;
 65:   Mat     A, B, C;

 67:   MatSchurComplementGetSubMatrices(pc->mat, &A, NULL, &B, &C, NULL);
 68:   KSPSolve(lsc->kspL, x, lsc->x1);
 69:   KSPCheckSolve(lsc->kspL, pc, lsc->x1);
 70:   MatMult(B, lsc->x1, lsc->x0);
 71:   if (lsc->scale) VecPointwiseMult(lsc->x0, lsc->x0, lsc->scale);
 72:   MatMult(A, lsc->x0, lsc->y0);
 73:   if (lsc->scale) VecPointwiseMult(lsc->y0, lsc->y0, lsc->scale);
 74:   MatMult(C, lsc->y0, lsc->x1);
 75:   KSPSolve(lsc->kspL, lsc->x1, y);
 76:   KSPCheckSolve(lsc->kspL, pc, y);
 77:   return 0;
 78: }

 80: static PetscErrorCode PCReset_LSC(PC pc)
 81: {
 82:   PC_LSC *lsc = (PC_LSC *)pc->data;

 84:   VecDestroy(&lsc->x0);
 85:   VecDestroy(&lsc->y0);
 86:   VecDestroy(&lsc->x1);
 87:   VecDestroy(&lsc->scale);
 88:   KSPDestroy(&lsc->kspL);
 89:   MatDestroy(&lsc->L);
 90:   return 0;
 91: }

 93: static PetscErrorCode PCDestroy_LSC(PC pc)
 94: {
 95:   PCReset_LSC(pc);
 96:   PetscFree(pc->data);
 97:   return 0;
 98: }

100: static PetscErrorCode PCSetFromOptions_LSC(PC pc, PetscOptionItems *PetscOptionsObject)
101: {
102:   PC_LSC *lsc = (PC_LSC *)pc->data;

104:   PetscOptionsHeadBegin(PetscOptionsObject, "LSC options");
105:   {
106:     PetscOptionsBool("-pc_lsc_scale_diag", "Use diagonal of velocity block (A) for scaling", "None", lsc->scalediag, &lsc->scalediag, NULL);
107:   }
108:   PetscOptionsHeadEnd();
109:   return 0;
110: }

112: static PetscErrorCode PCView_LSC(PC pc, PetscViewer viewer)
113: {
114:   PC_LSC   *jac = (PC_LSC *)pc->data;
115:   PetscBool iascii;

117:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
118:   if (iascii) {
119:     PetscViewerASCIIPushTab(viewer);
120:     if (jac->kspL) {
121:       KSPView(jac->kspL, viewer);
122:     } else {
123:       PetscViewerASCIIPrintf(viewer, "PCLSC KSP object not yet created, hence cannot display");
124:     }
125:     PetscViewerASCIIPopTab(viewer);
126:   }
127:   return 0;
128: }

130: /*MC
131:      PCLSC - Preconditioning for Schur complements, based on Least Squares Commutators

133:    Options Database Key:
134: .    -pc_lsc_scale_diag - Use the diagonal of A for scaling

136:    Level: intermediate

138:    Notes:
139:    This preconditioner will normally be used with `PCFIELDSPLIT` to precondition the Schur complement, but
140:    it can be used for any Schur complement system.  Consider the Schur complement

142: .vb
143:    S = A11 - A10 inv(A00) A01
144: .ve

146:    `PCLSC` currently doesn't do anything with A11, so let's assume it is 0.  The idea is that a good approximation to
147:    inv(S) is given by

149: .vb
150:    inv(A10 A01) A10 A00 A01 inv(A10 A01)
151: .ve

153:    The product A10 A01 can be computed for you, but you can provide it (this is
154:    usually more efficient anyway).  In the case of incompressible flow, A10 A01 is a Laplacian; call it L.  The current
155:    interface is to hang L and a preconditioning matrix Lp on the preconditioning matrix.

157:    If you had called `KSPSetOperators`(ksp,S,Sp), S should have type `MATSCHURCOMPLEMENT` and Sp can be any type you
158:    like (`PCLSC` doesn't use it directly) but should have matrices composed with it, under the names "LSC_L" and "LSC_Lp".
159:    For example, you might have setup code like this

161: .vb
162:    PetscObjectCompose((PetscObject)Sp,"LSC_L",(PetscObject)L);
163:    PetscObjectCompose((PetscObject)Sp,"LSC_Lp",(PetscObject)Lp);
164: .ve

166:    And then your Jacobian assembly would look like

168: .vb
169:    PetscObjectQuery((PetscObject)Sp,"LSC_L",(PetscObject*)&L);
170:    PetscObjectQuery((PetscObject)Sp,"LSC_Lp",(PetscObject*)&Lp);
171:    if (L) { assembly L }
172:    if (Lp) { assemble Lp }
173: .ve

175:    With this, you should be able to choose LSC preconditioning, using e.g. ML's algebraic multigrid to solve with L

177: .vb
178:    -fieldsplit_1_pc_type lsc -fieldsplit_1_lsc_pc_type ml
179: .ve

181:    Since we do not use the values in Sp, you can still put an assembled matrix there to use normal preconditioners.

183:    References:
184: +  * - Elman, Howle, Shadid, Shuttleworth, and Tuminaro, Block preconditioners based on approximate commutators, 2006.
185: -  * - Silvester, Elman, Kay, Wathen, Efficient preconditioning of the linearized Navier Stokes equations for incompressible flow, 2001.

187: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `Block_Preconditioners`, `PCFIELDSPLIT`,
188:           `PCFieldSplitGetSubKSP()`, `PCFieldSplitSetFields()`, `PCFieldSplitSetType()`, `PCFieldSplitSetIS()`, `PCFieldSplitSetSchurPre()`,
189:           `MatCreateSchurComplement()`, `MatCreateSchurComplement()`, `MatSchurComplementSetSubMatrices()`, `MatSchurComplementUpdateSubMatrices()`,
190:           `MatSchurComplementSetAinvType()`, `MatGetSchurComplement()`
191: M*/

193: PETSC_EXTERN PetscErrorCode PCCreate_LSC(PC pc)
194: {
195:   PC_LSC *lsc;

197:   PetscNew(&lsc);
198:   pc->data = (void *)lsc;

200:   pc->ops->apply           = PCApply_LSC;
201:   pc->ops->applytranspose  = NULL;
202:   pc->ops->setup           = PCSetUp_LSC;
203:   pc->ops->reset           = PCReset_LSC;
204:   pc->ops->destroy         = PCDestroy_LSC;
205:   pc->ops->setfromoptions  = PCSetFromOptions_LSC;
206:   pc->ops->view            = PCView_LSC;
207:   pc->ops->applyrichardson = NULL;
208:   return 0;
209: }