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: }