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:   PetscFunctionBegin;
 18:   if (lsc->allocated) PetscFunctionReturn(PETSC_SUCCESS);
 19:   PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &lsc->kspL));
 20:   PetscCall(KSPSetErrorIfNotConverged(lsc->kspL, pc->erroriffailure));
 21:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)lsc->kspL, (PetscObject)pc, 1));
 22:   PetscCall(KSPSetType(lsc->kspL, KSPPREONLY));
 23:   PetscCall(KSPSetOptionsPrefix(lsc->kspL, ((PetscObject)pc)->prefix));
 24:   PetscCall(KSPAppendOptionsPrefix(lsc->kspL, "lsc_"));
 25:   PetscCall(MatSchurComplementGetSubMatrices(pc->mat, &A, NULL, NULL, NULL, NULL));
 26:   PetscCall(MatCreateVecs(A, &lsc->x0, &lsc->y0));
 27:   PetscCall(MatCreateVecs(pc->pmat, &lsc->x1, NULL));
 28:   if (lsc->scalediag) PetscCall(VecDuplicate(lsc->x0, &lsc->scale));
 29:   lsc->allocated = PETSC_TRUE;
 30:   PetscFunctionReturn(PETSC_SUCCESS);
 31: }

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

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

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

 69:   PetscFunctionBegin;
 70:   PetscCall(MatSchurComplementGetSubMatrices(pc->mat, &A, NULL, &B, &C, NULL));
 71:   PetscCall(KSPSolve(lsc->kspL, x, lsc->x1));
 72:   PetscCall(KSPCheckSolve(lsc->kspL, pc, lsc->x1));
 73:   PetscCall(MatMult(B, lsc->x1, lsc->x0));
 74:   if (lsc->scale) PetscCall(VecPointwiseMult(lsc->x0, lsc->x0, lsc->scale));
 75:   PetscCall(MatMult(A, lsc->x0, lsc->y0));
 76:   if (lsc->scale) PetscCall(VecPointwiseMult(lsc->y0, lsc->y0, lsc->scale));
 77:   PetscCall(MatMult(C, lsc->y0, lsc->x1));
 78:   PetscCall(KSPSolve(lsc->kspL, lsc->x1, y));
 79:   PetscCall(KSPCheckSolve(lsc->kspL, pc, y));
 80:   PetscFunctionReturn(PETSC_SUCCESS);
 81: }

 83: static PetscErrorCode PCReset_LSC(PC pc)
 84: {
 85:   PC_LSC *lsc = (PC_LSC *)pc->data;

 87:   PetscFunctionBegin;
 88:   PetscCall(VecDestroy(&lsc->x0));
 89:   PetscCall(VecDestroy(&lsc->y0));
 90:   PetscCall(VecDestroy(&lsc->x1));
 91:   PetscCall(VecDestroy(&lsc->scale));
 92:   PetscCall(KSPDestroy(&lsc->kspL));
 93:   PetscCall(MatDestroy(&lsc->L));
 94:   PetscFunctionReturn(PETSC_SUCCESS);
 95: }

 97: static PetscErrorCode PCDestroy_LSC(PC pc)
 98: {
 99:   PetscFunctionBegin;
100:   PetscCall(PCReset_LSC(pc));
101:   PetscCall(PetscFree(pc->data));
102:   PetscFunctionReturn(PETSC_SUCCESS);
103: }

105: static PetscErrorCode PCSetFromOptions_LSC(PC pc, PetscOptionItems *PetscOptionsObject)
106: {
107:   PC_LSC *lsc = (PC_LSC *)pc->data;

109:   PetscFunctionBegin;
110:   PetscOptionsHeadBegin(PetscOptionsObject, "LSC options");
111:   {
112:     PetscCall(PetscOptionsBool("-pc_lsc_scale_diag", "Use diagonal of velocity block (A) for scaling", "None", lsc->scalediag, &lsc->scalediag, NULL));
113:   }
114:   PetscOptionsHeadEnd();
115:   PetscFunctionReturn(PETSC_SUCCESS);
116: }

118: static PetscErrorCode PCView_LSC(PC pc, PetscViewer viewer)
119: {
120:   PC_LSC   *jac = (PC_LSC *)pc->data;
121:   PetscBool iascii;

123:   PetscFunctionBegin;
124:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
125:   if (iascii) {
126:     PetscCall(PetscViewerASCIIPushTab(viewer));
127:     if (jac->kspL) {
128:       PetscCall(KSPView(jac->kspL, viewer));
129:     } else {
130:       PetscCall(PetscViewerASCIIPrintf(viewer, "PCLSC KSP object not yet created, hence cannot display"));
131:     }
132:     PetscCall(PetscViewerASCIIPopTab(viewer));
133:   }
134:   PetscFunctionReturn(PETSC_SUCCESS);
135: }

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

140:    Options Database Key:
141: .    -pc_lsc_scale_diag - Use the diagonal of A for scaling

143:    Level: intermediate

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

149: .vb
150:    S = A11 - A10 inv(A00) A01
151: .ve

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

156: .vb
157:    inv(A10 A01) A10 A00 A01 inv(A10 A01)
158: .ve

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

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

168: .vb
169:    PetscObjectCompose((PetscObject)Sp,"LSC_L",(PetscObject)L);
170:    PetscObjectCompose((PetscObject)Sp,"LSC_Lp",(PetscObject)Lp);
171: .ve

173:    And then your Jacobian assembly would look like

175: .vb
176:    PetscObjectQuery((PetscObject)Sp,"LSC_L",(PetscObject*)&L);
177:    PetscObjectQuery((PetscObject)Sp,"LSC_Lp",(PetscObject*)&Lp);
178:    if (L) { assembly L }
179:    if (Lp) { assemble Lp }
180: .ve

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

184: .vb
185:    -fieldsplit_1_pc_type lsc -fieldsplit_1_lsc_pc_type ml
186: .ve

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

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

194: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `Block_Preconditioners`, `PCFIELDSPLIT`,
195:           `PCFieldSplitGetSubKSP()`, `PCFieldSplitSetFields()`, `PCFieldSplitSetType()`, `PCFieldSplitSetIS()`, `PCFieldSplitSetSchurPre()`,
196:           `MatCreateSchurComplement()`, `MatCreateSchurComplement()`, `MatSchurComplementSetSubMatrices()`, `MatSchurComplementUpdateSubMatrices()`,
197:           `MatSchurComplementSetAinvType()`, `MatGetSchurComplement()`
198: M*/

200: PETSC_EXTERN PetscErrorCode PCCreate_LSC(PC pc)
201: {
202:   PC_LSC *lsc;

204:   PetscFunctionBegin;
205:   PetscCall(PetscNew(&lsc));
206:   pc->data = (void *)lsc;

208:   pc->ops->apply           = PCApply_LSC;
209:   pc->ops->applytranspose  = NULL;
210:   pc->ops->setup           = PCSetUp_LSC;
211:   pc->ops->reset           = PCReset_LSC;
212:   pc->ops->destroy         = PCDestroy_LSC;
213:   pc->ops->setfromoptions  = PCSetFromOptions_LSC;
214:   pc->ops->view            = PCView_LSC;
215:   pc->ops->applyrichardson = NULL;
216:   PetscFunctionReturn(PETSC_SUCCESS);
217: }