Actual source code: lsc.c

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

  3: typedef struct {
  4:   PetscBool allocated, commute, scalediag;
  5:   KSP       kspL, kspMass;
  6:   Vec       Avec0, Avec1, Svec0, scale;
  7:   Mat       L;
  8: } PC_LSC;

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

 15:   PetscFunctionBegin;
 16:   if (lsc->allocated) PetscFunctionReturn(PETSC_SUCCESS);
 17:   PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &lsc->kspL));
 18:   PetscCall(KSPSetNestLevel(lsc->kspL, pc->kspnestlevel));
 19:   PetscCall(KSPSetErrorIfNotConverged(lsc->kspL, pc->erroriffailure));
 20:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)lsc->kspL, (PetscObject)pc, 1));
 21:   PetscCall(KSPSetType(lsc->kspL, KSPPREONLY));
 22:   PetscCall(KSPSetOptionsPrefix(lsc->kspL, ((PetscObject)pc)->prefix));
 23:   PetscCall(KSPAppendOptionsPrefix(lsc->kspL, "lsc_"));
 24:   PetscCall(MatSchurComplementGetSubMatrices(pc->mat, &A, NULL, NULL, NULL, NULL));
 25:   PetscCall(MatCreateVecs(A, &lsc->Avec0, &lsc->Avec1));
 26:   PetscCall(MatCreateVecs(pc->pmat, &lsc->Svec0, NULL));
 27:   if (lsc->scalediag) PetscCall(VecDuplicate(lsc->Avec0, &lsc->scale));

 29:   if (lsc->commute) {
 30:     PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &lsc->kspMass));
 31:     PetscCall(KSPSetErrorIfNotConverged(lsc->kspMass, pc->erroriffailure));
 32:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)lsc->kspMass, (PetscObject)pc, 1));
 33:     PetscCall(KSPSetType(lsc->kspMass, KSPPREONLY));
 34:     PetscCall(KSPSetOptionsPrefix(lsc->kspMass, ((PetscObject)pc)->prefix));
 35:     PetscCall(KSPAppendOptionsPrefix(lsc->kspMass, "lsc_mass_"));
 36:   } else lsc->kspMass = NULL;

 38:   lsc->allocated = PETSC_TRUE;
 39:   PetscFunctionReturn(PETSC_SUCCESS);
 40: }

 42: static PetscErrorCode PCSetUp_LSC(PC pc)
 43: {
 44:   PC_LSC *lsc = (PC_LSC *)pc->data;
 45:   Mat     L, Lp, Qscale;

 47:   PetscFunctionBegin;
 48:   PetscCall(PCLSCAllocate_Private(pc));

 50:   /* Query for L operator */
 51:   PetscCall(PetscObjectQuery((PetscObject)pc->mat, "LSC_L", (PetscObject *)&L));
 52:   if (!L) PetscCall(PetscObjectQuery((PetscObject)pc->pmat, "LSC_L", (PetscObject *)&L));
 53:   PetscCall(PetscObjectQuery((PetscObject)pc->pmat, "LSC_Lp", (PetscObject *)&Lp));
 54:   if (!Lp) PetscCall(PetscObjectQuery((PetscObject)pc->mat, "LSC_Lp", (PetscObject *)&Lp));

 56:   /* Query for mass operator */
 57:   PetscCall(PetscObjectQuery((PetscObject)pc->pmat, "LSC_Qscale", (PetscObject *)&Qscale));
 58:   if (!Qscale) PetscCall(PetscObjectQuery((PetscObject)pc->mat, "LSC_Qscale", (PetscObject *)&Qscale));

 60:   if (lsc->commute) {
 61:     PetscCheck(L || Lp, PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "The user must provide an L operator for LSC preconditioning when commuting");
 62:     if (!L && Lp) L = Lp;
 63:     else if (L && !Lp) Lp = L;

 65:     PetscCheck(Qscale, PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "The user must provide a Qscale matrix for LSC preconditioning when commuting");
 66:   } else {
 67:     if (lsc->scale) {
 68:       if (!Qscale) PetscCall(MatSchurComplementGetSubMatrices(pc->mat, NULL, &Qscale, NULL, NULL, NULL));
 69:       PetscCall(MatGetDiagonal(Qscale, lsc->scale));
 70:       PetscCall(VecReciprocal(lsc->scale));
 71:     }
 72:     if (!L) {
 73:       Mat B, C;
 74:       PetscCall(MatSchurComplementGetSubMatrices(pc->mat, NULL, NULL, &B, &C, NULL));
 75:       if (lsc->scale) {
 76:         Mat CAdiaginv;
 77:         PetscCall(MatDuplicate(C, MAT_COPY_VALUES, &CAdiaginv));
 78:         PetscCall(MatDiagonalScale(CAdiaginv, NULL, lsc->scale));
 79:         if (!lsc->L) PetscCall(MatMatMult(CAdiaginv, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &lsc->L));
 80:         else PetscCall(MatMatMult(CAdiaginv, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &lsc->L));
 81:         PetscCall(MatDestroy(&CAdiaginv));
 82:       } else {
 83:         if (!lsc->L) {
 84:           PetscCall(MatProductCreate(C, B, NULL, &lsc->L));
 85:           PetscCall(MatProductSetType(lsc->L, MATPRODUCT_AB));
 86:           PetscCall(MatProductSetFromOptions(lsc->L));
 87:           PetscCall(MatProductSymbolic(lsc->L));
 88:         }
 89:         PetscCall(MatProductNumeric(lsc->L));
 90:       }
 91:       Lp = L = lsc->L;
 92:     }
 93:   }

 95:   PetscCall(KSPSetOperators(lsc->kspL, L, Lp));
 96:   PetscCall(KSPSetFromOptions(lsc->kspL));
 97:   if (lsc->commute) {
 98:     PetscCall(KSPSetOperators(lsc->kspMass, Qscale, Qscale));
 99:     PetscCall(KSPSetFromOptions(lsc->kspMass));
100:   }
101:   PetscFunctionReturn(PETSC_SUCCESS);
102: }

104: static PetscErrorCode PCApply_LSC(PC pc, Vec x, Vec y)
105: {
106:   PC_LSC *lsc = (PC_LSC *)pc->data;
107:   Mat     A, B, C;

109:   PetscFunctionBegin;
110:   PetscCall(MatSchurComplementGetSubMatrices(pc->mat, &A, NULL, &B, &C, NULL));
111:   if (lsc->commute) {
112:     PetscCall(KSPSolve(lsc->kspMass, x, lsc->Svec0));
113:     PetscCall(KSPCheckSolve(lsc->kspMass, pc, lsc->Svec0));
114:     PetscCall(MatMult(B, lsc->Svec0, lsc->Avec0));
115:     PetscCall(KSPSolve(lsc->kspL, lsc->Avec0, lsc->Avec1));
116:     PetscCall(KSPCheckSolve(lsc->kspL, pc, lsc->Avec1));
117:     PetscCall(MatMult(A, lsc->Avec1, lsc->Avec0));
118:     PetscCall(KSPSolve(lsc->kspL, lsc->Avec0, lsc->Avec1));
119:     PetscCall(KSPCheckSolve(lsc->kspL, pc, lsc->Avec1));
120:     PetscCall(MatMult(C, lsc->Avec1, lsc->Svec0));
121:     PetscCall(KSPSolve(lsc->kspMass, lsc->Svec0, y));
122:     PetscCall(KSPCheckSolve(lsc->kspMass, pc, y));
123:   } else {
124:     PetscCall(KSPSolve(lsc->kspL, x, lsc->Svec0));
125:     PetscCall(KSPCheckSolve(lsc->kspL, pc, lsc->Svec0));
126:     PetscCall(MatMult(B, lsc->Svec0, lsc->Avec0));
127:     if (lsc->scale) PetscCall(VecPointwiseMult(lsc->Avec0, lsc->Avec0, lsc->scale));
128:     PetscCall(MatMult(A, lsc->Avec0, lsc->Avec1));
129:     if (lsc->scale) PetscCall(VecPointwiseMult(lsc->Avec1, lsc->Avec1, lsc->scale));
130:     PetscCall(MatMult(C, lsc->Avec1, lsc->Svec0));
131:     PetscCall(KSPSolve(lsc->kspL, lsc->Svec0, y));
132:     PetscCall(KSPCheckSolve(lsc->kspL, pc, y));
133:   }
134:   PetscFunctionReturn(PETSC_SUCCESS);
135: }

137: static PetscErrorCode PCReset_LSC(PC pc)
138: {
139:   PC_LSC *lsc = (PC_LSC *)pc->data;

141:   PetscFunctionBegin;
142:   PetscCall(VecDestroy(&lsc->Avec0));
143:   PetscCall(VecDestroy(&lsc->Avec1));
144:   PetscCall(VecDestroy(&lsc->Svec0));
145:   PetscCall(KSPDestroy(&lsc->kspL));
146:   if (lsc->commute) PetscCall(KSPDestroy(&lsc->kspMass));
147:   if (lsc->L) PetscCall(MatDestroy(&lsc->L));
148:   if (lsc->scale) PetscCall(VecDestroy(&lsc->scale));
149:   PetscFunctionReturn(PETSC_SUCCESS);
150: }

152: static PetscErrorCode PCDestroy_LSC(PC pc)
153: {
154:   PetscFunctionBegin;
155:   PetscCall(PCReset_LSC(pc));
156:   PetscCall(PetscFree(pc->data));
157:   PetscFunctionReturn(PETSC_SUCCESS);
158: }

160: static PetscErrorCode PCSetFromOptions_LSC(PC pc, PetscOptionItems *PetscOptionsObject)
161: {
162:   PC_LSC *lsc = (PC_LSC *)pc->data;

164:   PetscFunctionBegin;
165:   PetscOptionsHeadBegin(PetscOptionsObject, "LSC options");
166:   {
167:     PetscCall(PetscOptionsBool("-pc_lsc_commute", "Whether to commute the LSC preconditioner in the style of Olshanskii", "None", lsc->commute, &lsc->commute, NULL));
168:     PetscCall(PetscOptionsBool("-pc_lsc_scale_diag", "Whether to scale BBt products. Will use the inverse of the diagonal of Qscale or A if the former is not provided.", "None", lsc->scalediag, &lsc->scalediag, NULL));
169:     PetscCheck(!lsc->scalediag || !lsc->commute, PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "Diagonal-based scaling is not used when doing a commuted LSC. Either do not ask for diagonal-based scaling or use non-commuted LSC in the original style of Elman");
170:   }
171:   PetscOptionsHeadEnd();
172:   PetscFunctionReturn(PETSC_SUCCESS);
173: }

175: static PetscErrorCode PCView_LSC(PC pc, PetscViewer viewer)
176: {
177:   PC_LSC   *jac = (PC_LSC *)pc->data;
178:   PetscBool iascii;

180:   PetscFunctionBegin;
181:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
182:   if (iascii) {
183:     PetscCall(PetscViewerASCIIPushTab(viewer));
184:     if (jac->kspL) {
185:       PetscCall(KSPView(jac->kspL, viewer));
186:     } else {
187:       PetscCall(PetscViewerASCIIPrintf(viewer, "PCLSC KSP object not yet created, hence cannot display"));
188:     }
189:     if (jac->commute) {
190:       if (jac->kspMass) {
191:         PetscCall(KSPView(jac->kspMass, viewer));
192:       } else {
193:         PetscCall(PetscViewerASCIIPrintf(viewer, "PCLSC Mass KSP object not yet created, hence cannot display"));
194:       }
195:     }
196:     PetscCall(PetscViewerASCIIPopTab(viewer));
197:   }
198:   PetscFunctionReturn(PETSC_SUCCESS);
199: }

201: /*MC
202:    PCLSC - Preconditioning for Schur complements, based on Least Squares Commutators {cite}`elmanhowleshadidshuttleworthtuminaro2006` {cite}`silvester2001efficient`

204:    Options Database Key:
205: +    -pc_lsc_commute    - Whether to commute the LSC preconditioner in the style of Olshanskii
206: -    -pc_lsc_scale_diag - Whether to scale $BB^T$ products. Will use the inverse of the diagonal of Qscale or A if the former is not provided

208:    Level: intermediate

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

214:    $$
215:    S = A11 - A10 A00^{-1} A01
216:    $$

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

221:    $$
222:    (A10 A01)^{-1} A10 A00 A01 (A10 A01)^{-1}
223:    $$

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

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

233: .vb
234:    PetscObjectCompose((PetscObject)Sp,"LSC_L",(PetscObject)L);
235:    PetscObjectCompose((PetscObject)Sp,"LSC_Lp",(PetscObject)Lp);
236: .ve

238:    And then your Jacobian assembly would look like

240: .vb
241:    PetscObjectQuery((PetscObject)Sp,"LSC_L",(PetscObject*)&L);
242:    PetscObjectQuery((PetscObject)Sp,"LSC_Lp",(PetscObject*)&Lp);
243:    if (L) { assembly L }
244:    if (Lp) { assemble Lp }
245: .ve

247:    With this, you should be able to choose LSC preconditioning, using e.g. the `PCML` algebraic multigrid to solve with L
248: .vb
249:    -fieldsplit_1_pc_type lsc -fieldsplit_1_lsc_pc_type ml
250: .ve

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

254: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `Block_Preconditioners`, `PCFIELDSPLIT`,
255:           `PCFieldSplitGetSubKSP()`, `PCFieldSplitSetFields()`, `PCFieldSplitSetType()`, `PCFieldSplitSetIS()`, `PCFieldSplitSetSchurPre()`,
256:           `MatCreateSchurComplement()`, `MatCreateSchurComplement()`, `MatSchurComplementSetSubMatrices()`, `MatSchurComplementUpdateSubMatrices()`,
257:           `MatSchurComplementSetAinvType()`, `MatGetSchurComplement()`
258: M*/

260: PETSC_EXTERN PetscErrorCode PCCreate_LSC(PC pc)
261: {
262:   PC_LSC *lsc;

264:   PetscFunctionBegin;
265:   PetscCall(PetscNew(&lsc));
266:   pc->data = (void *)lsc;

268:   pc->ops->apply           = PCApply_LSC;
269:   pc->ops->applytranspose  = NULL;
270:   pc->ops->setup           = PCSetUp_LSC;
271:   pc->ops->reset           = PCReset_LSC;
272:   pc->ops->destroy         = PCDestroy_LSC;
273:   pc->ops->setfromoptions  = PCSetFromOptions_LSC;
274:   pc->ops->view            = PCView_LSC;
275:   pc->ops->applyrichardson = NULL;
276:   PetscFunctionReturn(PETSC_SUCCESS);
277: }